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

Free Energy, Rates, and Mechanism of Transmembrane Dimerization in Lipid Bilayers from Dynamically Unbiased Molecular Dynamics Simulations

The assembly of proteins in membranes plays a key role in many crucial cellular pathways. Despite their importance, characterizing transmembrane assembly remains challenging for experiments and simulations. Equilibrium molecular dynamics simulations do not cover the time scales required to sample the typical transmembrane assembly. Hence, most studies rely on enhanced sampling schemes that steer the dynamics of transmembrane proteins along a collective variable that should encode all slow degrees of freedom. However, given the complexity of the condensed-phase lipid environment, this is far from trivial, with the consequence that free energy profiles of dimerization can be poorly converged. Here, we introduce an alternative approach, which relies only on simulating short, dynamically unbiased paths, avoiding using collective variables or biasing forces. By merging all paths, we obtain free energy profiles, rates, and mechanisms of transmembrane dimerization with the same set of simulations. We showcase our algorithm by sampling the spontaneous association and dissociation of a transmembrane protein in a lipid bilayer, the popular coarse-grained Martini force field. Our algorithm represents a promising way to investigate assembly processes in biologically relevant membranes, overcoming some of the challenges of conventional methods.

Characterization of hyperactive mutations in the renal potassium channel ROMK uncovers unique effects on channel biogenesis and ion conductance

Hypertension affects one billion people worldwide and is the most common risk factor for cardiovascular disease, yet a comprehensive picture of its underlying genetic factors is incomplete. Amongst regulators of blood pressure is the renal outer medullary potassium (ROMK) channel. While select ROMK mutants are prone to premature degradation and lead to disease, heterozygous carriers of some of these same alleles are protected from hypertension. Therefore, we hypothesized that gain-of-function (GoF) ROMK variants which increase potassium flux may predispose people to hypertension. To begin to test this hypothesis, we employed genetic screens and a candidate-based approach to identify six GoF variants in yeast. Subsequent functional assays in higher cells revealed two variant classes. The first group exhibited greater stability in the endoplasmic reticulum, enhanced channel assembly, and/or increased protein at the cell surface. The second group of variants resided in the PIP2-binding pocket, and computational modeling coupled with patch-clamp studies demonstrated lower free energy for channel opening and slowed current rundown, consistent with an acquired PIP2-activated state. Together, these findings advance our understanding of ROMK structure-function, suggest the existence of hyperactive ROMK alleles in humans, and establish a system to facilitate the development of ROMK-targeted antihypertensives.

Machine Learning Derived Collective Variables for the Study of Protein Homodimerization in Membrane

The accurate calculation of equilibrium constants for protein-protein association is of fundamental importance to quantitative biology and remains an outstanding challenge for computational biophysics. Traditionally, equilibrium constants have been computed from one-dimensional free energy surfaces derived from sampling along a single collective variable. Importantly, recent advances in enhanced sampling methodology have facilitated the characterization of multidimensional free energy landscapes, often exposing multiple thermodynamically important minima missed by more restrictive sampling methods. A key to the effectiveness of this multidimensional sampling approach is the identification of collective variables that effectively define the configurational space of dissociated and associated states. Here we present the application of two machine learning methods for the unbiased determination of collective variables for enhanced sampling and analysis of protein-protein association. Our results both validate prior work, based on intuition derived collective variables, and demonstrate the effectiveness of the machine learning methods for the identification of collective variables for association reactions in complex biomolecular systems.

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.

Preferential electrostatic interactions of phosphatidic acid with arginines

Phosphatidic acid (PA) is an anionic lipid that preferentially interacts with proteins in a diverse set of cellular processes such as transport, apoptosis, and neurotransmission. One such interaction is that of the PA lipids with the proteins of voltage-sensitive ion channels. In comparison to several other similarly charged anionic lipids, PA lipids exhibit much stronger interactions. Intrigued and motivated by this finding, we sought out to gain deeper understanding into the electrostatic interactions of anionic lipids with charged proteins. Using the voltage sensor domain (VSD) of the KvAP channel as a model system, we performed long-timescale atomistic simulations to analyze the interactions of POPA, POPG, and POPI lipids with arginines (ARGs). Our simulations reveal two mechanisms. First, POPA is able to interact not only with surface ARGs but is able to snorkel and interact with a buried arginine. POPG and POPI lipids on the other hand show weak interactions even with both the surface and buried ARGs. Second, deprotonated POPA with -2 charge is able to break the salt-bridge connection between VSD protein segments and establish its own electrostatic bond with the ARG. Based on these findings, we propose a headgroup size hypothesis for preferential solvation of proteins by charged lipids. These findings may be valuable in understanding how PA lipids could be modulating kinematics of transmembrane proteins in cellular membranes.

Mechanism of hydrophobic gating in the acetylcholine receptor channel pore

Neuromuscular acetylcholine receptors (AChRs) are hetero-pentameric, ligand-gated ion channels. The binding of the neurotransmitter acetylcholine (ACh) to two target sites promotes a global conformational change of the receptor that opens the channel and allows ion conduction through the channel pore. Here, by measuring free-energy changes from single-channel current recordings and using molecular dynamics simulations, we elucidate how a constricted hydrophobic region acts as a "gate" to regulate the channel opening in the pore of AChRs. Mutations of gate residues, including those implicated in congenital myasthenia syndrome, lower the permeation barrier of the channel substantially and increase the unliganded gating equilibrium constant (constitutive channel openings). Correlations between hydrophobicity and the observed free-energy changes, supported by calculations of water densities in the wild-type versus mutant channel pores, provide evidence for hydrophobic wetting-dewetting transition at the gate. The analysis of a coupled interaction network provides insight into the molecular mechanism of closed- versus open-state conformational changes at the gate. Studies of the transition state by "phi"(φ)-value analysis indicate that agonist binding serves to stabilize both the transition and the open state. Intersubunit interaction energy measurements and molecular dynamics simulations suggest that channel opening involves tilting of the pore-lining M2 helices, asymmetric outward rotation of amino acid side chains, and wetting transition of the gate region that lowers the barrier to ion permeation and stabilizes the channel open conformation. Our work provides new insight into the hydrophobic gate opening and shows why the gate mutations result in constitutive AChR channel activity.

Bilayer tension-induced clustering of the UPR sensor IRE1

The endoplasmic reticulum acts as a protein quality control center where a range of chaperones and foldases facilitates protein folding. IRE1 is a sensory transmembrane protein that transduces signals of proteotoxic stress by forming clusters and activating a cellular program called the unfolded protein response (UPR). Recently, membrane thickness variation due to membrane compositional changes have been shown to drive IRE1 cluster formation, activating the UPR even in the absence of proteotoxic stress. Here, we demonstrate a direct relationship between bilayer tension and UPR activation based on IRE1 dimer stability. The stability of the IRE1 dimer in a (50%DOPC-50%POPC) membrane at different applied bilayer tensions was analyzed via molecular dynamics simulations. The potential of mean force for IRE1 dimerization predicts a higher concentration of IRE1 dimers for both tensed and compressed ER membranes. This study shows that IRE1 may be a mechanosensitive membrane protein and establishes a direct biophysical relationship between bilayer tension and UPR activation.

Effect of AQP4 and its palmitoylation on the permeability of exogenous reactive oxygen species: Insights from computational study

Aquaporin 4 (AQP4) facilitates the transport of reactive oxygen species (ROS). Both cancer cells and the ionizing radiation microenvironment can induce posttranslational modifications (PTMs) in AQP4, which may affect its permeability to ROS. Because this ROS diffusion process is rapid, microscopic, and instantaneous within and outside cells, conventional experimental methods are inadequate for elucidating the molecular mechanisms involved. In this study, computational methods were employed to investigate the permeability of exogenous ROS mediated by radiation in AQP4 at a molecular scale. We constructed a simulation system incorporating AQP4 and AQP4-Cysp13 in a complex lipid environment with ROS. Long-timescale molecular dynamics simulations were conducted to assess the structural stability of both AQP4 and AQP4-Cysp13. Free energy calculations were utilized to determine the ROS transport capability of the two AQP4 proteins. Computational electrophysiology and channel structural analysis quantitatively evaluated changes in ROS transport capacity under various radiation-induced transmembrane voltage microenvironments. Our findings demonstrate the distinct transport capabilities of AQP4 channels for water molecules and various types of ROS and reveal a decrease in transport efficiency when AQP4 undergoes palmitoylation modification. In addition, we have simulated the radiation-induced alteration of cell membrane voltage, which significantly affected the ROS transport capacity. We propose that this research will enhance the understanding of the molecular mechanisms governing the transport of exogenous ROS by AQP4 and elucidate the influence of palmitoylation on ROS transport. This study will also help clarify how different structural features of AQP4 affect the transport of exogenous ROS mediated by radiotherapy, thereby providing a theoretical molecular basis for the development of new treatment strategies that combine with radiotherapy.

A Rigorous Framework for Calculating Protein-Protein Binding Affinities in Membranes

Calculating the binding free energy of integral transmembrane (TM) proteins is crucial for understanding the mechanisms by which they recognize one another and reversibly associate. The glycophorin A (GpA) homodimer, composed of two α-helical segments, has long served as a model system for studying TM protein reversible association. The present work establishes a methodological framework for calculating the binding affinity of the GpA homodimer in the heterogeneous environment of a membrane. Our investigation carefully considered a variety of protocols, including the appropriate choice of the force field, rigorous standardization reflecting the experimental conditions, sampling algorithm, anisotropic environment, and collective variables, to accurately describe GpA dimerization via molecular dynamics-based approaches. Specifically, two strategies were explored: (i) an unrestrained potential mean force (PMF) calculation, which merely enhances sampling along the separation of the two binding partners without any restraint, and (ii) a so-called "geometrical route", whereby the α-helices are progressively separated with imposed restraints on their orientational, positional, and conformational degrees of freedom to accelerate convergence. Our simulations reveal that the simplified, unrestrained PMF approach is inadequate for the description of GpA dimerization. Instead, the geometrical route, tailored specifically to GpA in a membrane environment, yields excellent agreement with experimental data within a reasonable computational time. A dimerization free energy of -10.7 kcal/mol is obtained, in fairly good agreement with available experimental data. The geometrical route further helps elucidate how environmental forces drive association before helical interactions stabilize it. Our simulations also brought to light a distinct, long-lived spatial arrangement that potentially serves as an intermediate state during dimer formation. The methodological advances in the generalized geometrical route provide a powerful tool for accurate and efficient binding-affinity calculations of intricate TM protein complexes in inhomogeneous environments.

The role of structural heterogeneity in the homodimerization of transmembrane proteins

The equilibrium association of transmembrane proteins plays a fundamental role in membrane protein function and cellular signaling. While the study of the equilibrium binding of single pass transmembrane proteins has received significant attention in experiment and simulation, the accurate assessment of equilibrium association constants remains a challenge to experiment and simulation. In experiment, there remain wide variations in association constants derived from experimental studies of the most widely studied transmembrane proteins. In simulation, state-of-the art methods have failed to adequately sample the thermodynamically relevant structures of the dimer state ensembles using coarse-grained models. In addition, all-atom force fields often fail to accurately assess the relative free energies of the dimer and monomer states. Given the importance of this fundamental biophysical process, it is essential to address these shortcomings. In this work, we establish an effective computational protocol for the calculation of equilibrium association constants for transmembrane homodimer formation. A set of transmembrane protein homodimers, used in the parameterization of the MARTINI v3 force field, are simulated using metadynamics, based on three collective variables. The method is found to be accurate and computationally efficient, providing a standard to be used in the future simulation studies using coarse-grained or all-atom models.

Umbrella Sampling Simulations of Cardiac Thin Filament Reveal Thermodynamic Consequences of Troponin I Inhibitory Peptide Mutations

The cardiac thin filament comprises F-actin, tropomyosin, and troponin (cTn). cTn is composed of three subunits: troponin C (cTnC), troponin I (cTnI), and troponin T (cTnT). To computationally study the effect of the thin filament on cTn activation events, we employed targeted molecular dynamics followed by umbrella sampling using a model of the thin filament to measure the thermodynamics of cTn transition events. Our simulations revealed that the thin filament causes an increase in the free energy required to open the cTnC hydrophobic patch and causes a more favorable interaction between this region and the cTnI switch peptide. Mutations to the cTn complex can lead to cardiomyopathy, a collection of diseases that present clinically with symptoms of hypertrophy or dilation of the cardiac muscle, leading to impairment of the heart's ability to function normally and ultimately myocardial infarction or heart failure. Upon introduction of cardiomyopathic mutations to R145 of cTnI, we observed a general decrease in the free energy of opening the cTnC hydrophobic patch, which is on par with previous experimental results. These mutations also exhibited a decrease in electrostatic interactions between cTnI-R145 and actin-E334. After introduction of a small molecule to the wild-type cTnI-actin interface to intentionally disrupt intersubunit contacts, we successfully observed similar thermodynamic consequences and disruptions to the same protein-protein contacts as observed with the cardiomyopathic mutations. Computational studies utilizing the cTn complex in isolation would have been unable to observe these effects, highlighting the importance of using a more physiologically relevant thin-filament model to investigate the global consequences of cardiomyopathic mutations to the cTn complex.

Structural insights into Nirmatrelvir (PF-07321332)-3C-like SARS-CoV-2 protease complexation: a ligand Gaussian accelerated molecular dynamics study

Coronavirus 3C-like protease (3CLpro) is found in SARS-CoV-2 virus, which causes COVID-19. 3CLpro controls virus replication and is a major target for target-based antiviral discovery. As reported by Pfizer, Nirmatrelvir (PF-07321332) is a competitive protein inhibitor and a clinical candidate for orally delivered medication. However, the binding mechanisms between Nirmatrelvir and 3CLpro complex structures remain unknown. This study incorporated ligand Gaussian accelerated molecular dynamics, the one-dimensional and two-dimensional potential of mean force, normal molecular dynamics, and Kramers' rate theory to determine the binding and dissociation rate constants (k_off and k_on) associated with the binding of the 3CLpro protein to the Nirmatrelvir inhibitor. The proposed approach addresses the challenges in designing small-molecule antiviral drugs.

The guidance and adhesion protein FLRT2 dimerizes in cis via dual small-X3-small transmembrane motifs

Fibronectin Leucine-rich Repeat Transmembrane (FLRT 1-3) proteins are a family of broadly expressed single-spanning transmembrane receptors that play key roles in development. Their extracellular domains mediate homotypic cell-cell adhesion and heterotypic protein interactions with other receptors to regulate cell adhesion and guidance. These in trans FLRT interactions determine the formation of signaling complexes of varying complexity and function. Whether FLRTs also interact at the surface of the same cell, in cis, remains unknown. Here, molecular dynamics simulations reveal two dimerization motifs in the FLRT2 transmembrane helix. Single particle tracking experiments show that these Small-X₃-Small motifs synergize with a third dimerization motif encoded in the extracellular domain to permit the cis association and co-diffusion patterns of FLRT2 receptors on cells. These results may point to a competitive switching mechanism between in cis and in trans interactions, which suggests that homotypic FLRT interaction mirrors the functionalities of classic adhesion molecules.

Molecular Mechanism of CD44 Homodimerization Modulated by Palmitoylation and Membrane Environments

The homodimerization of CD44 plays a key role in an intercellular-to-extracellular signal transduction and tumor progression. Acylated modification and specific membrane environments have been reported to mediate translocation and oligomerization of CD44, however, the underlying molecular mechanism remains elusive. In this study, extensive molecular dynamics simulations are performed to characterize the dimerization of palmitoylated CD44 variants in different bilayer environments. CD44 forms homodimer depending on the cysteines on the juxta-membrane domains, and the dimerization efficiency and packing configurations are defected by their palmitoylated modifications. In the phase-segregated (raft included) membrane, homodimerization of the palmitoylated CD44 is hardly observed, whereas PIP2 addition compensates to realize dimerization. However, the structure of CD44 homodimer formed in the phase-segregated bilayer turns susceptive and PIP2 addition allows for an extensive conformation of the cytoplasmic domain, a proposal prerequisite to access the cytoskeleton linker proteins. The results unravel a delicate competitive relationship between PIP2, lipid raft and palmitoylation in mediating protein homodimerization, which helps to clarify the dynamic dimer conformations and involved cellular signaling of the CD44 likewise proteins.

On Computing Equilibrium Binding Constants for Protein-Protein Association in Membranes

Protein association in lipid membranes is fundamental to membrane protein function and of great biomedical relevance. All-atom and coarse-grained models have been extensively used to understand the protein-protein interactions in the membrane and to compute equilibrium association constants. However, slow translational and rotational diffusion of protein in membrane presents challenges to the effective sampling of conformations defining the ensembles of free and bound states contributing to the association equilibrium and the free energy of dimerization. We revisit the homodimerization equilibrium of the TM region of glycophorin A. Conformational sampling is performed using umbrella sampling along previously proposed one-dimensional collective variables and compared with sampling over a two-dimensional collective variable space using the MARTINI v2.2 force field. We demonstrate that the one-dimensional collective variables suffer from restricted sampling of the native homodimer conformations leading to a biased free energy landscape. Conversely, simulations along the two-dimensional collective variable effectively characterize the thermodynamically relevant native and non-native interactions contributing to the association equilibrium. These results demonstrate the challenges associated with accurately characterizing binding equilibria when multiple poses contribute to the bound state ensemble.

Effects of Cholesterol and PIP2 on Interactions between Glycophorin A and Band 3 in Lipid Bilayers

In the erythrocyte membrane, the interactions between Glycophorin A (GPA) and Band 3 are associated strongly with the biological function of the membrane and several blood disorders. In this work, using coarse-grained molecular dynamics simulation, we systematically investigate the effects of cholesterol and phosphatidylinositol-4,5-bisphosphate (PIP2) on the interactions of GPA with Band 3 in the model erythrocyte membranes. We examine the dynamics of the interactions of GPA with Band 3 in different lipid bilayers on the microsecond time scale and calculate the binding free energy between GPA and Band 3. The results indicate that cholesterols thermodynamically favor the binding of GPA to Band 3 by increasing the thickness of the lipid bilayer and by producing an effective attraction between the proteins due to the depletion effect. Cholesterols also slow the kinetics of the binding of GPA to Band 3 by reducing the lateral mobility of the lipids and proteins and may influence the binding sites between the proteins. The anionic PIP2 lipids prefer binding to the surface of the proteins through electrostatic attraction between the PIP2 headgroup and the positively charged residues on the protein surface. Ions in the solvent facilitate the PIP2 aggregation which promotes the binding of GPA to Band 3.

Lipid-Mediated Association of the Slg1 Transmembrane Domains in Yeast Plasma Membranes

Clustering of transmembrane proteins underlies a multitude of fundamental biological processes at the plasma membrane (PM) such as receptor activation, lateral domain formation, and mechanotransduction. The self-association of the respective transmembrane domains (TMDs) has also been suggested to be responsible for the micron-scaled patterns seen for integral membrane proteins in the budding yeast PM. However, the underlying interplay between the local lipid composition and the TMD identity is still not mechanistically understood. In this work, we combined coarse-grained molecular dynamics simulations of simplified bilayer systems with high-resolution live-cell microscopy to analyze the distribution of a representative helical yeast TMD from the PM sensor Slg1 within different lipid environments. In our simulations, we specifically evaluated the effects of acyl chain saturation and anionic lipid head groups on the association of two TMDs. We found that weak lipid-protein interactions significantly affect the configuration of TMD dimers and the free energy of association. Increased amounts of unsaturated phospholipids (PLs) strongly reduced the helix-helix interaction, while the presence of anionic phosphatidylserine (PS) hardly affected the dimer formation. We could experimentally confirm this surprising lack of effect of PS using the network factor, a mesoscopic measure of PM pattern formation in yeast cells. Simulations also showed that the formation of TMD dimers in turn increased the order parameter of the surrounding lipids and induced long-range perturbations in lipid organization. In summary, our results shed new light on the mechanisms of lipid-mediated dimerization of TMDs in complex lipid mixtures.

Probing the Formation, Structure and Free Energy Relationships of M Protein Dimers of SARS-CoV-2

The M protein of the novel coronavirus 2019 (SARS-CoV-2) is the major structural component of the viral envelope and is also the minimum requirement for virus particle budding. M proteins generally exist as dimers. In virus assembly, they are the main driving force for envelope formation through lateral interactions and interactions with other viral structural proteins that play a central role. We built 100 candidate models and finally analyzed the six most convincing structural features of the SARS-CoV-2 M protein dimer based on long-timescale molecular dynamics (MD) simulations, multiple free energy analyses (potential mean force (PMF) and molecular mechanics Poisson-Boltzmann surface area (MMPBSA)) and principal component analysis (PCA) to obtain the most reasonable structure. The dimer stability was found to depend on the Leu-Ile zipper motif and aromatic amino acids in the transmembrane domain (TMD). Furthermore, the C-terminal domain (CTD) effects were relatively small. These results highlight a model in which there is sufficient binding affinity between the TMDs of M proteins to form dimers through the residues at the interface of the three transmembrane helices (TMHs). This study aims to help find more effective inhibitors of SARS-CoV-2 M dimers and to develop vaccines based on structural information.

Substrate Access Mechanism in a Novel Membrane-Bound Phospholipase A of Pseudomonas aeruginosa Concordant with Specificity and Regioselectivity

PlaF is a cytoplasmic membrane-bound phospholipase A₁ from Pseudomonas aeruginosa that alters the membrane glycerophospholipid (GPL) composition and fosters the virulence of this human pathogen. PlaF activity is regulated by a dimer-to-monomer transition followed by tilting of the monomer in the membrane. However, how substrates reach the active site and how the characteristics of the active site tunnels determine the activity, specificity, and regioselectivity of PlaF for natural GPL substrates have remained elusive. Here, we combined unbiased and biased all-atom molecular dynamics (MD) simulations and configurational free-energy computations to identify access pathways of GPL substrates to the catalytic center of PlaF. Our results map out a distinct tunnel through which substrates access the catalytic center. PlaF variants with bulky tryptophan residues in this tunnel revealed decreased catalysis rates due to tunnel blockage. The MD simulations suggest that GPLs preferably enter the active site with the sn-1 acyl chain first, which agrees with the experimentally demonstrated PLA₁ activity of PlaF. We propose that the acyl chain-length specificity of PlaF is determined by the structural features of the access tunnel, which results in favorable free energy of binding of medium-chain GPLs. The suggested egress route conveys fatty acid (FA) products to the dimerization interface and, thus, contributes to understanding the product feedback regulation of PlaF by FA-triggered dimerization. These findings open up opportunities for developing potential PlaF inhibitors, which may act as antibiotics against P. aeruginosa.

The impact of phosphatidylserine exposure on cancer cell membranes on the activity of the anticancer peptide HB43

HB43 (FAKLLAKLAKKLL) is a synthetic peptide active against cell lines derived from breast, colon, melanoma, lung, prostate, and cervical cancers. Despite its remarkable spectrum of activity, the mechanism of action at the molecular level has never been investigated, preventing further optimization of its selectivity. The alternation of charged and hydrophobic residues suggests amphipathicity, but the formation of alpha-helical structure seems discouraged by its short length and the large number of positively charged residues. Using different biophysical and in silico approaches we show that HB43 is completely unstructured in solution but assumes alpha-helical conformation in the presence of DPC micelles and liposomes exposing phosphatidylserine (PS) used as mimics of cancer cell membranes. Membrane permeabilization assays demonstrate that the interaction leads to the preferential destabilization of PS-containing vesicles with respect to PC-containing ones, here used as noncancerous cell mimics. ssNMR reveals that HB43 is able to fluidify the internal structure of cancer-cell mimicking liposomes while MD simulations show its internalization in such bilayers. This is achieved by the formation of specific interactions between the lysine side chains and the carboxylate group of phosphatidylserine and/or the phosphate oxygen atoms of targeted phospholipids, which could catalyze the formation of the alpha helix required for internalization. With the aim of better understanding the peptide biocompatibility and the additional antibacterial activity, the interaction with noncancerous cell mimicking liposomes exposing phosphatidylcholine (PC) and bacterial mimicking bilayers exposing phosphatidylglycerol (PG) is also described.

From Bench to Biomolecular Simulation: Phospholipid Modulation of Potassium Channels

Potassium (K+) ion channels are crucial in numerous cellular processes as they hyperpolarise a cell through K+ conductance, returning a cell to its resting potential. K+ channel mutations result in multiple clinical complications such as arrhythmia, neonatal diabetes and migraines. Since 1995, the regulation of K+ channels by phospholipids has been heavily studied using a range of interdisciplinary methods such as cellular electrophysiology, structural biology and computational modelling. As a result, K+ channels are model proteins for the analysis of protein-lipid interactions. In this review, we will focus on the roles of lipids in the regulation of K+ channels, and how atomic-level structures, along with experimental techniques and molecular simulations, have helped guide our understanding of the importance of phospholipid interactions.

A double bilayer to study the nonequilibrium environmental response of GIRK2 in complex states

G protein-gated inwardly rectifying potassium (GIRK) channels play essential roles in electrical signaling in neurons and muscle cells. Nonequilibrium environments provide crucial driving forces behind many cellular events. Here, we apply the antiparallel alignment double bilayer model to study GIRK2 in response to the time-dependent membrane potential. Using molecular dynamics and umbrella sampling, we examined the time-dependent environmental impact on the ion conduction, energy basis, and primary motions of GIRK2 in different complex states with phosphatidylinositol-4,5-bisphosphate (PIP2) and G-protein βγ subunits (Gβγ). The antiparallel alignment double bilayer model enables us to study the transport performance in inward and outward K+ and mixed K+ and Na+. We obtained the recoverable discharge process of GIRK2 complexed with both PIP2 and Gβγ, compared with occasional conduction under PIP2-only regulation. Calculations of potential of mean force suggest different regulation by the helix bundle crossing (HBC) gate and G-loop gate regarding different complex states and under a membrane potential. In a nonequilibrium environment, distinct functional rocking motions of GIRK2 were identified under strengthened correlations between the transmembrane helices and downstream cytoplasmic domains with binding of PIP2, cations, and Gβγ. The findings suggest the potential domain motions and dynamics associated with a nonequilibrium environment and highlight the application of the antiparallel alignment double bilayer model to investigate factors in an asymmetric environment.

Building intuition for binding free energy calculations: Bound state definition, restraints, and symmetry

The theory behind computation of absolute binding free energies using explicit-solvent molecular simulations is well-established, yet somewhat complex, with counter-intuitive aspects. This leads to frequent frustration, common misconceptions, and sometimes erroneous numerical treatment. To improve this, we present the main practically relevant segments of the theory with constant reference to physical intuition. We pinpoint the role of the implicit or explicit definition of the bound state (or the binding site) to make a robust link between an experimental measurement and a computational result. We clarify the role of symmetry and discuss cases where symmetry number corrections have been misinterpreted. In particular, we argue that symmetry corrections as classically presented are a source of confusion and could be advantageously replaced by restraint free energy contributions. We establish that contrary to a common intuition, partial or missing sampling of some modes of symmetric bound states does not affect the calculated decoupling free energies. Finally, we review these questions and pitfalls in the context of a few common practical situations: binding to a symmetric receptor (equivalent binding sites), binding of a symmetric ligand (equivalent poses), and formation of a symmetric complex, in the case of homodimerization.

Multiple dimerizing motifs at different locations modulate the dimerization of the syndecan transmembrane domains

Syndecans (SDCs) are a family of four members of integral membrane proteins, which play important roles in cell-cell interactions. Dimerization/oligomerization generated by transmembrane domains (TMDs) appears to crucially regulate several functional behaviors of all syndecan members. The different levels of protein-protein interactions mediated by Syndecan TMDs may lead to a rather complicated function of Syndecans. The molecular mechanism of the different dimerization tendencies in each type of SDCs remains unclear. Here, the self-assembly process of syndecan TMD homodimers and heterodimers was studied in molecular details by molecular dynamics simulations. Our computational results showed that the SDC2 forms the most stable homodimer, which is consistent with previous experimental results. Detailed analysis suggests that instead of the conserved dimerizing motif G8XXXG12 in all four SDCs involved in homo- and hetero-dimerization of SDCs. The different locations of GXXXA motif affect the stability of SDC dimers. In addition, we found that A3XXXA7 can stabilize the dimerization, making the dimer of SDC2 the most stable among these SDC dimers. Our results shed light on the complex effect of multiple dimerizing motifs on the dimerization of transmembrane domains.

Molecular simulation studies of the interactions between the human/pangolin/cat/bat ACE2 and the receptor binding domain of the SARS-CoV-2 spike protein

The recent outbreak of SARS-CoV-2 has had a profound effect on the world. Similar to that in SARS-CoV, the entry receptor of SARS-CoV-2 is ACE2. The binding of SARS-CoV-2 spike protein to ACE2 is the critical to the virus infection. Recently multiple species (human, Chinese chrysanthemum, Malay pangolin and cat) have been reported to be susceptible to the virus infection. However, the binding capacity and the detailed binding mechanism of SARS-CoV-2 spike protein to ACE2 of these species remains unexplored. Herein free energy calculations with MM-GBSA and Potential of Mean Forces together reveal that the Human-SARS-CoV-2 has a higher stability tendency than Human-SARS-CoV. Meanwhile, we uncover that SARS-CoV-2 has an enhanced ability to bind with the ACE2 in humans, pangolins and cats compared to that in bats. Analysis of key residues with energy decomposition and residue contact maps reveal several important consensus sites in ACE2s among the studied species, and determined the more favorable specified residues among the different types of amino acids. These results provide important implications for understanding SARS-CoV-2 host range which will make it possible to control the spread of the virus and use of animal models, targeted drug screening and vaccine candidates against SARS-CoV-2.

Coarse-Grained Modeling of Multiple Pathways in Conformational Transitions of Multi-Domain Proteins

Large-scale conformational transitions in multi-domain proteins are often essential for their functions. To investigate the transitions, it is necessary to explore multiple potential pathways, which involve different intermediate structures. Here, we present a multi-basin (MB) coarse-grained (CG) structure-based Go̅ model for describing transitions in proteins with more than two moving domains. This model is an extension of our dual-basin Go̅ model in which system-dependent parameters are determined systematically using the multistate Bennett acceptance ratio method. In the MB Go̅ model for multi-domain proteins, we assume that intermediate structures may have partial inter-domain native contacts. This approach allows us to search multiple transition pathways that involve distinct intermediate structures using the CG molecular dynamics (MD) simulations. We apply this scheme to an enzyme, adenylate kinase (AdK), which has three major domains and can move along two different pathways. Using the optimized mixing parameters for each pathway, AdK shows frequent transitions between the Open, Closed, and the intermediate basins and samples a wide variety of conformations within each basin. The explored multiple transition pathways could be compared with experimental data and examined in more detail by atomistic MD simulations.

Antimicrobial Bombinin-like Peptide 3 Selectively Recognizes and Inserts into Bacterial Biomimetic Bilayers in Multiple Steps

Bombinins are a wide family of antimicrobial peptides from Xenopus skin. By sequence clustering, we highlighted at least three families named A, B, and H, which might exert antibacterial activity by different modes of action. In this work, we study bombinin-like peptide 3 (BLP-3) as a nonhemolytic representative of the quite unexplored class A due to its appealing activity toward WHO-priority-list bacteria such as Neisseria, Pseudomonas aeruginosa, and Staphylococcus aureus. A marked preference for cardiolipin and phosphatidylglycerol head groups, typically found in bacteria, is proven with biomimetic membranes studied by liquid and solid NMR and MD simulations. BLP-3 gets structured upon interaction and penetrates deeply into the bilayer in two steps involving a superficial insertion of key side chains and subsequent internalization. All along the pathway, a fundamental role is played by lysine residues in the conserved region 11-19, which act in synergy with other key residues.

Molecular dynamics simulations of bacterial outer membrane lipid extraction: Adequate sampling?

The outer membrane of Gram-negative bacteria is almost exclusively composed of lipopolysaccharide in its outer leaflet, whereas the inner leaflet contains a mixture of phospholipids. Lipopolysaccharide diffuses at least an order of magnitude slower than phospholipids, which can cause issues for molecular dynamics simulations in terms of adequate sampling. Here, we test a number of simulation protocols for their ability to achieve convergence with reasonable computational effort using the MARTINI coarse-grained force-field. This is tested in the context both of potential of mean force (PMF) calculations for lipid extraction from membranes and of lateral mixing within the membrane phase. We find that decoupling the cations that cross-link the lipopolysaccharide headgroups from the extracted lipid during PMF calculations is the best approach to achieve convergence comparable to that for phospholipid extraction. We also show that lateral lipopolysaccharide mixing/sorting is very slow and not readily addressable even with Hamiltonian replica exchange. We discuss why more sorting may be unrealistic for the short (microseconds) timescales we simulate and provide an outlook for future studies of lipopolysaccharide-containing membranes.

Interaction energy between neuronal cell receptors and peptoid ligands

Peptoids as an extracellular matrix (ECM) material is gaining importance in in vitro neuronal cell culture studies due to their biocompatibility, self-assembling structure, and stability. Mechanotransduction between a neuronal cell and an ECM is mediated by neuronal cell receptors such as integrin and neural cellular adhesion molecule. In this study, using molecular dynamics, we investigate the interaction energies between peptoid and neuronal cell receptors, and also study the effect of peptoid bundle size. We investigate the interaction surface between peptoid bundles and neuronal cell receptors, integrin and neural cellular adhesion molecule, using the solvent accessible surface area method to find the influence of hydrophobic and hydrophilic residues of the peptoid chain. We find the free energy landscape using the umbrella sampling method and then evaluate the potential mean force (PMF) and unbinding force during the dissociation between peptoid bundles and neuronal cell receptors. We find that the peptoid bundles have a higher affinity for the neuronal cell receptors, however increasing the size of peptoid bundles increases the affinity for integrin and neural cell adhesion molecule. PMF data for peptoid and neuronal cell receptor dissociation indicates that binding force increases as the size of the peptoid bundle increases. The higher binding strength during peptoid and neuronal cell receptors are due to the hydrophobic residue cluster area in the binding region. These findings will provide a better insight into using peptoid as an ECM.

Comparison of umbrella sampling and steered molecular dynamics methods for computing free energy profiles of aromatic substrates through phospholipid bilayers

Understanding the permeation of molecules through lipid membranes is fundamental for predicting the cellular uptake of solutes and drug delivery mechanisms. In molecular simulations, the usual approach is to compute the free energy (FE) profile of a molecule across a model lipid bilayer, which can then be used to estimate the permeability of the molecule. Umbrella Sampling (US), which involves carrying out a series of biased simulations along a defined reaction coordinate (usually the bilayer normal direction), is a popular method for the computation of such FE profiles. However, US can be challenging to implement because the results are dependent on the strength of the biasing potential and the spacing of windows along the reaction coordinate, which, in practice, are usually optimized by an inefficient trial and error approach. The Steered Molecular Dynamics implementation of the Jarzynski Equality (JE-SMD) has been identified as an alternative to equilibrium sampling methods for measuring the FE change across a reaction coordinate. In the JE-SMD approach, equilibrium FE values are evaluated from the average of rapid non-equilibrium trajectories, thus avoiding the practical issues that come with US. Here, we use three different corrections of the JE-SMD method to calculate the FE change for the translocation of two aromatic substrates, phenylalanine and toluene, across a lipid bilayer and compare the accuracy and computational efficiency of these approaches to the results obtained using US. We show evidence that when computing the free energy profile, the JE-SMD approach suffers from insufficient sampling convergence of the bilayer environment and is dependent on the characteristic of the aromatic substrate itself. We deduce that, despite its drawbacks, US remains the more viable approach of the two for computing the FE profile.

Molecular Basis of the Anticancer and Antibacterial Properties of CecropinXJ Peptide: An In Silico Study

Esophageal cancer is an aggressive lethal malignancy causing thousands of deaths every year. While current treatments have poor outcomes, cecropinXJ (CXJ) is one of the very few peptides with demonstrated in vivo activity. The great interest in CXJ stems from its low toxicity and additional activity against most ESKAPE bacteria and fungi. Here, we present the first study of its mechanism of action based on molecular dynamics (MD) simulations and sequence-property alignment. Although unstructured in solution, predictions highlight the presence of two helices separated by a flexible hinge containing P24 and stabilized by the interaction of W2 with target biomembranes: an amphipathic helix-I and a poorly structured helix-II. Both MD and sequence-property alignment point to the important role of helix I in both the activity and the interaction with biomembranes. MD reveals that CXJ interacts mainly with phosphatidylserine (PS) but also with phosphatidylethanolamine (PE) headgroups, both found in the outer leaflet of cancer cells, while salt bridges with phosphate moieties are prevalent in bacterial biomimetic membranes composed of PE, phosphatidylglycerol (PG) and cardiolipin (CL). The antibacterial activity of CXJ might also explain its interaction with mitochondria, whose phospholipid composition recalls that of bacteria and its capability to induce apoptosis in cancer cells.

Guardians of the Cell: State-of-the-Art of Membrane Proteins from a Computational Point-of-View

Membrane proteins (MPs) encompass a large family of proteins with distinct cellular functions, and although representing over 50% of existing pharmaceutical drug targets, their structural and functional information is still very scarce. Over the last years, in silico analysis and algorithm development were essential to characterize MPs and overcome some limitations of experimental approaches. The optimization and improvement of these methods remain an ongoing process, with key advances in MPs' structure, folding, and interface prediction being continuously tackled. Herein, we discuss the latest trends in computational methods toward a deeper understanding of the atomistic and mechanistic details of MPs.

Molecular Dynamics Simulations in Drug Discovery and Pharmaceutical Development

Molecular dynamics (MD) simulations have become increasingly useful in the modern drug development process. In this review, we give a broad overview of the current application possibilities of MD in drug discovery and pharmaceutical development. Starting from the target validation step of the drug development process, we give several examples of how MD studies can give important insights into the dynamics and function of identified drug targets such as sirtuins, RAS proteins, or intrinsically disordered proteins. The role of MD in antibody design is also reviewed. In the lead discovery and lead optimization phases, MD facilitates the evaluation of the binding energetics and kinetics of the ligand-receptor interactions, therefore guiding the choice of the best candidate molecules for further development. The importance of considering the biological lipid bilayer environment in the MD simulations of membrane proteins is also discussed, using G-protein coupled receptors and ion channels as well as the drug-metabolizing cytochrome P450 enzymes as relevant examples. Lastly, we discuss the emerging role of MD simulations in facilitating the pharmaceutical formulation development of drugs and candidate drugs. Specifically, we look at how MD can be used in studying the crystalline and amorphous solids, the stability of amorphous drug or drug-polymer formulations, and drug solubility. Moreover, since nanoparticle drug formulations are of great interest in the field of drug delivery research, different applications of nano-particle simulations are also briefly summarized using multiple recent studies as examples. In the future, the role of MD simulations in facilitating the drug development process is likely to grow substantially with the increasing computer power and advancements in the development of force fields and enhanced MD methodologies.

Combining Virtual Reality Visualization with Ensemble Molecular Dynamics to Study Complex Protein Conformational Changes

Molecular dynamics (MD) simulations are increasingly used to elucidate relationships between protein structure, dynamics, and their biological function. Currently, it is extremely challenging to perform MD simulations of large-scale structural rearrangements in proteins that occur on millisecond timescales or beyond, as this requires very significant computational resources, or the use of cumbersome "collective variable" enhanced sampling protocols. Here, we describe a framework that combines ensemble MD simulations and virtual reality visualization (eMD-VR) to enable users to interactively generate realistic descriptions of large amplitude, millisecond timescale protein conformational changes in proteins. Detailed tests demonstrate that eMD-VR substantially decreases the computational cost of folding simulations of a WW domain, without the need to define collective variables a priori. We further show that eMD-VR generated pathways can be combined with Markov state models to describe the thermodynamics and kinetics of large-scale loop motions in the enzyme cyclophilin A. Our results suggest eMD-VR is a powerful tool for exploring protein energy landscapes in bioengineering efforts.

Antimicrobial Peptide K11 Selectively Recognizes Bacterial Biomimetic Membranes and Acts by Twisting Their Bilayers

K11 is a synthetic peptide originating from the introduction of a lysine residue in position 11 within the sequence of a rationally designed antibacterial scaffold. Despite its remarkable antibacterial properties towards many ESKAPE bacteria and its optimal therapeutic index (320), a detailed description of its mechanism of action is missing. As most antimicrobial peptides act by destabilizing the membranes of the target organisms, we investigated the interaction of K11 with biomimetic membranes of various phospholipid compositions by liquid and solid-state NMR. Our data show that K11 can selectively destabilize bacterial biomimetic membranes and torque the surface of their bilayers. The same is observed for membranes containing other negatively charged phospholipids which might suggest additional biological activities. Molecular dynamic simulations reveal that K11 can penetrate the membrane in four steps: after binding to phosphate groups by means of the lysine residue at the N-terminus (anchoring), three couples of lysine residues act subsequently to exert a torque in the membrane (twisting) which allows the insertion of aromatic side chains at both termini (insertion) eventually leading to the flip of the amphipathic helix inside the bilayer core (helix flip and internalization).

Evaluating inositol phospholipid interactions with inward rectifier potassium channels and characterising their role in disease

Membrane proteins are frequently modulated by specific protein-lipid interactions. The activation of human inward rectifying potassium (hKir) channels by phosphoinositides (PI) has been well characterised. Here, we apply a coarse-grained molecular dynamics free-energy perturbation (CG-FEP) protocol to capture the energetics of binding of PI lipids to hKir channels. By using either a single- or multi-step approach, we establish a consistent value for the binding of PIP2 to hKir channels, relative to the binding of the bulk phosphatidylcholine phospholipid. Furthermore, by perturbing amino acid side chains on hKir6.2, we show that the neonatal diabetes mutation E179K increases PIP2 affinity, while the congenital hyperinsulinism mutation K67N results in a reduced affinity. We show good agreement with electrophysiological data where E179K exhibits a reduction in neomycin sensitivity, implying that PIP2 binds more tightly E179K channels. This illustrates the application of CG-FEP to compare affinities between lipid species, and for annotating amino acid residues.

Assessing the Perturbing Effects of Drugs on Lipid Bilayers Using Gramicidin Channel-Based In Silico and In Vitro Assays

Partitioning of bioactive molecules, including drugs, into cell membranes may produce indiscriminate changes in membrane protein function. As a guide to safe drug development, it therefore becomes important to be able to predict the bilayer-perturbing potency of hydrophobic/amphiphilic drugs candidates. Toward this end, we exploited gramicidin channels as molecular force probes and developed in silico and in vitro assays to measure drugs' bilayer-modifying potency. We examined eight drug-like molecules that were found to enhance or suppress gramicidin channel function in a thick 1,2-dierucoyl-sn-glycero-3-phosphocholine (DC22:1PC) but not in thin 1,2-dioleoyl-sn-glycero-3-phosphocholine (DC18:1PC) lipid bilayer. The mechanism underlying this difference was attributable to the changes in gramicidin dimerization free energy by drug-induced perturbations of lipid bilayer physical properties and bilayer-gramicidin interactions. The combined in silico and in vitro approaches, which allow for predicting the perturbing effects of drug candidates on membrane protein function, have implications for preclinical drug safety assessment.

Water Dynamics Around Proteins: T- and R-States of Hemoglobin and Melittin

The water dynamics, as characterized by the local hydrophobicity (LH), is investigated for tetrameric hemoglobin (Hb) and dimeric melittin. For the T₀ to R₀ transition in Hb, it is found that LH provides additional molecular-level insight into the Perutz mechanism, i.e., the breaking and formation of salt bridges at the α₁/β₂ and α₂/β₁ interface is accompanied by changes in LH. For Hb in cubic water boxes with 90 and 120 Å edge length it is observed that following a decrease in LH as a consequence of reduced water density or change of water orientation at the protein/water interface the α/β interfaces are destabilized; this is a hallmark of the Perutz stereochemical model for the T to R transition in Hb. The present work thus provides a dynamical view of the classical structural model relevant to the molecular foundations of Hb function. For dimeric melittin, earlier results by Cheng and Rossky [ Nature 1998, 392, 696-699] are confirmed and interpreted on the basis of LH from simulations in which the protein structure is frozen. For the flexible melittin dimer, the changes in the local hydration can be as much as 30% greater than for the rigid dimer, reflecting the fact that protein and water dynamics are coupled.

Atomistic mechanism of transmembrane helix association

Transmembrane helix association is a fundamental step in the folding of helical membrane proteins. The prototypical example of this association is formation of the glycophorin dimer. While its structure and stability have been well-characterized experimentally, the detailed assembly mechanism is harder to obtain. Here, we use all-atom simulations within phospholipid membrane to study glycophorin association. We find that initial association results in the formation of a non-native intermediate, separated by a significant free energy barrier from the dimer with a native binding interface. We have used transition-path sampling to determine the association mechanism. We find that the mechanism of the initial bimolecular association to form the intermediate state can be mediated by many possible contacts, but seems to be particularly favoured by formation of non-native contacts between the C-termini of the two helices. On the other hand, the contacts which are key to determining progression from the intermediate to the native state are those which define the native binding interface, reminiscent of the role played by native contacts in determining folding of globular proteins. As a check on the simulations, we have computed association and dissociation rates from the transition-path sampling. We obtain results in reasonable accord with available experimental data, after correcting for differences in native state stability. Our results yield an atomistic description of the mechanism for a simple prototype of helical membrane protein folding.

Many important cellular functions are performed by membrane proteins, and in particular by association of proteins via transmembrane helices. However, the mechanism of how the helices associate has been challenging to study, by either experiment or simulation. Here, we use advanced molecular simulation methods to overcome the slow time scales involved in helix association and dissociation and obtain a view of the association mechanism in atomic detail. We show that association occurs via an initially non-native dimer, before proceeding to the native state, and we validate our results by comparison to available experimental kinetic data. Our methods will also aid in the study of the assembly mechanism of larger transmembrane proteins via molecular simulation.

Characterizing Membrane Association and Periplasmic Transfer of Bacterial Lipoproteins through Molecular Dynamics Simulations

Escherichia coli lipoprotein precursors at the inner membrane undergo three maturation stages before transport by the Lol system to the outer membrane. Here, we develop a pipeline to simulate the membrane association of bacterial lipoproteins in their four maturation states. This has enabled us to model and simulate 81 of the predicted 114 E. coli lipoproteins and reveal their interactions with the host lipid membrane. As part of this set we characterize the membrane contacts of LolB, the lipoprotein involved in periplasmic translocation. We also consider the means and bioenergetics for lipoprotein localization. Our calculations uncover a preference for LolB over LolA and therefore indicate how a lipoprotein may be favorably transferred from the inner to outer membrane. Finally, we reveal that LolC has a role in membrane destabilization, thereby promoting lipoprotein transfer to LolA.

Inserting Small Molecules across Membrane Mixtures: Insight from the Potential of Mean Force

Small solutes have been shown to alter the lateral organization of cell membranes and reconstituted phospholipid bilayers; however, the mechanisms by which these changes happen are still largely unknown. Traditionally, both experiment and simulation studies have been restricted to testing only a few compounds at a time, failing to identify general molecular descriptors or chemical properties that would allow extrapolating beyond the subset of considered solutes. In this work, we probe the competing energetics of inserting a solute in different membrane environments by means of the potential of mean force. We show that these calculations can be used as a computationally efficient proxy to establish whether a solute will stabilize or destabilize domain phase separation. Combined with umbrella-sampling simulations and coarse-grained molecular dynamics simulations, we are able to screen solutes across a wide range of chemistries and polarities. Our results indicate that for the system under consideration, preferential partitioning and therefore effectiveness in altering membrane phase separation are strictly linked to the location of insertion in the bilayer (i.e., midplane or interface). Our approach represents a fast and simple tool for obtaining structural and thermodynamic insight into the partitioning of small molecules between lipid domains and its relation to phase separation, ultimately providing a platform for identifying the key determinants of this process.

The energetics of protein-lipid interactions as viewed by molecular simulations

Membranes are formed from a bilayer containing diverse lipid species with which membrane proteins interact. Integral, membrane proteins are embedded in this bilayer, where they interact with lipids from their surroundings, whilst peripheral membrane proteins bind to lipids at the surface of membranes. Lipid interactions can influence the function of membrane proteins, either directly or allosterically. Both experimental (structural) and computational approaches can reveal lipid binding sites on membrane proteins. It is, therefore, important to understand the free energies of these interactions. This affords a more complete view of the engagement of a particular protein with the biological membrane surrounding it. Here, we describe many computational approaches currently in use for this purpose, including recent advances using both free energy and unbiased simulation methods. In particular, we focus on interactions of integral membrane proteins with cholesterol, and with anionic lipids such as phosphatidylinositol 4,5-bis-phosphate and cardiolipin. Peripheral membrane proteins are exemplified via interactions of PH domains with phosphoinositide-containing membranes. We summarise the current state of the field and provide an outlook on likely future directions of investigation.

Enhanced sampling and free energy calculations for protein simulations

Molecular dynamics simulation is a powerful computational technique to study biomolecular systems, which complements experiments by providing insights into the structural dynamics relevant to biological functions at atomic scale. It can also be used to calculate the free energy landscapes of the conformational transitions to better understand the functions of the biomolecules. However, the sampling of biomolecular configurations is limited by the free energy barriers that need to be overcome, leading to considerable gaps between the timescales reached by MD simulation and those governing biological processes. To address this issue, many enhanced sampling methodologies have been developed to increase the sampling efficiency of molecular dynamics simulations and free energy calculations. Usually, enhanced sampling algorithms can be classified into methods based on collective variables (CV-based) and approaches which do not require predefined CVs (CV-free). In this chapter, the theoretical basis of free energy estimation is briefly reviewed first, followed by the reviews of the most common CV-based and CV-free methods including the presentation of some examples and recent developments. Finally, the combination of different enhanced sampling methods is discussed.

Multiple lipid binding sites determine the affinity of PH domains for phosphoinositide-containing membranes

Association of peripheral proteins with lipid bilayers regulates membrane signaling and dynamics. Pleckstrin homology (PH) domains bind to phosphatidylinositol phosphate (PIP) molecules in membranes. The effects of local PIP enrichment on the interaction of PH domains with membranes is unclear. Molecular dynamics simulations allow estimation of the binding energy of GRP1 PH domain to PIP3-containing membranes. The free energy of interaction of the PH domain with more than two PIP3 molecules is comparable to experimental values, suggesting that PH domain binding involves local clustering of PIP molecules within membranes. We describe a mechanism of PH binding proceeding via an encounter state to two bound states which differ in the orientation of the protein relative to the membrane, these orientations depending on the local PIP concentration. These results suggest that nanoscale clustering of PIP molecules can control the strength and orientation of PH domain interaction in a concentration-dependent manner.

Lipid Interactions of a Ciliary Membrane TRP Channel: Simulation and Structural Studies of Polycystin-2

Polycystin-2 (PC2) is a transient receptor potential (TRP) channel present in ciliary membranes of the kidney. PC2 shares a transmembrane fold with other TRP channels, in addition to an extracellular domain found in TRPP and TRPML channels. Using molecular dynamics (MD) simulations and cryoelectron microscopy we identify and characterize PIP₂ and cholesterol interactions with PC2. PC2 is revealed to have a PIP binding site close to the equivalent vanilloid/lipid binding site in the TRPV1 channel. A 3.0-Å structure reveals a binding site for cholesterol on PC2. Cholesterol interactions with the channel at this site are characterized by MD simulations. The two classes of lipid binding sites are compared with sites observed in other TRPs and in Kv channels. These findings suggest PC2, in common with other ion channels, may be modulated by both PIPs and cholesterol, and position PC2 within an emerging model of the roles of lipids in the regulation and organization of ciliary membranes.

•

Lipid interactions of PC2 channels have been explored by MD simulation and cryo-EM

•

PIP₂ binds to a site corresponding to the vanilloid/lipid binding site of TRPV1

•

Cholesterol binds between the S3 and S4 helices and S6 of the adjacent subunit

•

PC2, in common with other channels, may be modulated by PIPs and cholesterol

Building a macro-mixing dual-basin Gō model using the Multistate Bennett Acceptance Ratio

The dual-basin Gō-model is a structural-based coarsegrained model for simulating a conformational transition between two known structures of a protein. Two parameters are required to produce a dual-basin potential mixed using two single-basin potentials, although the determination of mixing parameters is usually not straightforward. Here, we have developed an efficient scheme to determine the mixing parameters using the Multistate Bennett Acceptance Ratio (MBAR) method after short simulations with a set of parameters. In the scheme, MBAR allows us to predict observables at various unsimulated conditions, which are useful to improve the mixing parameters in the next round of iterative simulations. The number of iterations that are necessary for obtaining the converged mixing parameters are significantly reduced in the scheme. We applied the scheme to two proteins, the glutamine binding protein and the ribose binding protein, for showing the effectiveness in the parameter determination. After obtaining the converged parameters, both proteins show frequent conformational transitions between open and closed states, providing the theoretical basis to investigate structure-dynamics-function relationships of the proteins.