One Descriptor to Fold Them All: Harnessing Intuition and Machine Learning to Identify Transferable Lasso Peptide Reaction Coordinates

Identifying optimal reaction coordinates for complex conformational changes and protein folding remains an outstanding challenge. This study combines collective variable (CV) discovery based on chemical intuition and machine learning with enhanced sampling to converge the folding free energy landscape of lasso peptides, a unique class of natural products with knot-like tertiary structures. This knotted scaffold imparts remarkable stability, making lasso peptides resistant to proteolytic degradation, thermal denaturation, and extreme pH conditions. Although their direct synthesis would enable therapeutic design, it has not yet been possible due to the improbable occurrence of spontaneous lasso folding. Thus, simulations characterizing the folding propensity are needed to identify strategies for increasing access to the lasso architecture by stabilizing the pre-lasso ensemble before isopeptide bond formation. Herein, harmonic linear discriminant analysis (HLDA) is combined with metadynamics-enhanced sampling to discover CVs capable of distinguishing the pre-lasso fold and converging the folding propensity. Intuitive CVs are compared to iterative rounds of HLDA to identify CVs that not only accomplish these goals for one lasso peptide but also seem to be transferable to others, establishing a protocol for the identification of folding reaction coordinates for lasso peptides.

Lasso Peptides: Exploring the Folding Landscape of Nature's Smallest Interlocked Motifs

Lasso peptides make up a class of natural products characterized by a threaded structure. Given their small size and stability, chemical synthesis would offer tremendous potential for the development of novel therapeutics. However, the accessibility of the pre-folded lasso architecture has limited this advance. To better understand the folding process de novo, simulations are used herein to characterize the folding propensity of microcin J25 (MccJ25), a lasso peptide known for its antimicrobial properties. New algorithms are developed to unambiguously distinguish threaded from nonthreaded precursors and determine handedness, a key feature in natural lasso peptides. We find that MccJ25 indeed forms right-handed pre-lassos, in contrast to past predictions but consistent with all natural lasso peptides. Additionally, the native pre-lasso structure is shown to be metastable prior to ring formation but to readily transition to entropically favored unfolded and nonthreaded structures, suggesting that de novo lasso folding is rare. However, by altering the ring forming residues and appending thiol and thioester functionalities, we are able to increase the stability of pre-lasso conformations. Furthermore, conditions leading to protonation of a histidine imidazole side chain further stabilize the modified pre-lasso ensemble. This work highlights the use of computational methods to characterize lasso folding and demonstrates that de novo access to lasso structures can be facilitated by optimizing sequence, unnatural modifications, and reaction conditions like pH.

tICA-Metadynamics for Identifying Slow Dynamics in Membrane Permeation

Molecular simulations are commonly used to understand the mechanism of membrane permeation of small molecules, particularly for biomedical and pharmaceutical applications. However, despite significant advances in computing power and algorithms, calculating an accurate permeation free energy profile remains elusive for many drug molecules because it can require identifying the rate-limiting degrees of freedom (i.e., appropriate reaction coordinates). To resolve this issue, researchers have developed machine learning approaches to identify slow system dynamics. In this work, we apply time-lagged independent component analysis (tICA), an unsupervised dimensionality reduction algorithm, to molecular dynamics simulations with well-tempered metadynamics to find the slowest collective degrees of freedom of the permeation process of trimethoprim through a multicomponent membrane. We show that tICA-metadynamics yields translational and orientational collective variables (CVs) that increase convergence efficiency ∼1.5 times. However, crossing the periodic boundary is shown to introduce artifacts in the translational CV that can be corrected by taking absolute values of molecular features. Additionally, we find that the convergence of the tICA CVs is reached with approximately five membrane crossings and that data reweighting is required to avoid deviations in the translational CV.

tICA-Metadynamics for Identifying Slow Dynamics in Membrane Permeation

Molecular simulations are commonly used to understand the mechanism of membrane permeation of small molecules, particularly for biomedical and pharmaceutical applications. However, despite significant advances in computing power and algorithms, calculating an accurate permeation free energy profile remains elusive for many drug molecules because it can require identifying the rate-limiting degrees of freedom (i.e., appropriate reaction coordinates). To resolve this issue, researchers have developed machine learning approaches to identify slow system dynamics. In this work, we apply time-lagged independent component analysis (tICA), an unsupervised dimensionality reduction algorithm, to molecular dynamics simulations with well-tempered metadynamics to find the slowest collective degrees of freedom of the permeation process of trimethoprim through a multicomponent membrane. We show that tICA-metadynamics yields translational and orientational collective variables (CVs) that increase convergence efficiency ∼1.5 times. However, crossing the periodic boundary is shown to introduce artefacts in the translational CV that can be corrected by taking absolute values of molecular features. Additionally, we find that the convergence of the tICA CVs is reached with approximately five membrane crossings, and that data reweighting is required to avoid deviations in the translational CV.

Mycolactone A vs. B: Multiscale Simulations Reveal the Roles of Localization and Association in Isomer-Specific Toxicity

Mycolactone is an exotoxin produced by Mycobacterium ulcerans that causes the neglected tropical skin disease Buruli ulcer. This toxin inhibits the Sec61 translocon in the endoplasmic reticulum (ER), preventing the host cell from producing several secretory and transmembrane proteins, resulting in cytotoxic and immunomodulatory effects. Interestingly, only one of the two dominant isoforms of mycolactone is cytotoxic. Here, we investigate the origin of this specificity by performing extensive molecular dynamics (MD) simulations with enhanced free energy sampling to query the association trends of the two isoforms with both the Sec61 translocon, using two distinct cryo-electron microscopy (cryo-EM) models as references, and the ER membrane, which serves as a toxin reservoir prior to association. Our results suggest that mycolactone B (the cytotoxic isoform) has a stronger association with the ER membrane than mycolactone A due to more favorable interactions with membrane lipids and water molecules. This could increase the reservoir of toxin proximal to the Sec61 translocon. In one model of Sec61 inhibited by mycolactone, we find that isomer B interacts more closely with residues thought to play a key role in signal peptide recognition and, thus, are essential for subsequent protein translocation. In the other model, we find that isomer B interacts more closely with the lumenal and lateral gates of the translocon, the dynamics of which are essential for protein translocation. These interactions induce a more closed conformation, which has been suggested to block signal peptide insertion and subsequent protein translocation. Collectively, these findings suggest that isomer B's unique cytotoxicity is a consequence of both increased localization to the ER membrane and channel-locking association with the Sec61 translocon, facets that could be targeted in the development of Buruli Ulcer diagnostics and Sec61-targeted therapeutics.

Mycolactone A vs. B: Does localization or association explain isomer-specific toxicity?

Mycolactone is an exotoxin produced by Mycobacterium ulcerans that causes the neglected tropical skin disease Buruli ulcer. This toxin inhibits the Sec61 translocon in the endoplasmic reticulum (ER), preventing the host cell from producing many secretory and transmembrane proteins, resulting in cytotoxic and immunomodulatory effects. Interestingly, only one of the two dominant isoforms of mycolactone is cytotoxic. Here, we investigate the origin of this specificity by performing extensive molecular dynamics (MD) simulations with enhanced free energy sampling to query the association trends of the two isoforms with both the Sec61 translocon and the ER membrane, which serves as a toxin reservoir prior to association. Our results suggest that mycolactone B (the cytotoxic isoform) has a stronger association with the ER membrane than mycolactone A due to more favorable interactions with membrane lipids and water molecules. This could increase the reservoir of toxin proximal to the Sec61 translocon. Isomer B also interacts more closely with the lumenal and lateral gates of the translocon, the dynamics of which are essential for protein translocation. These interactions induce a more closed conformation, which has been suggested to block signal peptide insertion and subsequent protein translocation. Collectively, these findings suggest that isomer B's unique cytotoxicity is a consequence of both increased localization to the ER membrane and channel-locking association with the Sec61 translocon, facets that could be targeted in the development of Buruli Ulcer diagnostics and Sec61-targeted therapeutics.

Can membrane composition traffic toxins? Mycolactone and preferential membrane interactions

Mycolactone is a cytotoxic and immunosuppressive macrolide produced by Mycobacterium ulcerans and the sole causative agent of the neglected tropical skin disease Buruli ulcer. The toxin acts by invading host cells and interacting with intracellular targets to disrupt multiple fundamental cellular processes. Mycolactone's amphiphilic nature enables strong interactions with lipophilic environments, including cellular membranes; however, the specificity of these interactions and the role of membranes in the toxin's pathogenicity remain unknown. It is likely that preferential interactions with lipophilic carriers play a key role in the toxin's distribution in the host, which, if understood, could provide insights to aid in the development of needed diagnostics for Buruli ulcer disease. In this work, molecular dynamics simulations were combined with enhanced free-energy sampling to characterize mycolactone's association with and permeation through models of the mammalian endoplasmic reticulum (ER) and plasma membranes (PMs). We find that increased order in the PMs not only leads to a different permeation mechanism compared with that in the ER membrane but also an energetic driving force for ER localization. Increased hydration, membrane deformation, and preferential interactions with unsaturated lipid tails stabilize the toxin in the ER membrane, while disruption of lipid packing is a destabilizing force in the PMs.

Capturing the Liquid-Crystalline Phase Transformation: Implications for Protein Targeting to Sterol Ester-Rich Lipid Droplets

Lipid droplets are essential organelles that store and traffic neutral lipids. The phospholipid monolayer surrounding their neutral lipid core engages with a highly dynamic proteome that changes according to cellular and metabolic conditions. Recent work has demonstrated that when the abundance of sterol esters increases above a critical concentration, such as under conditions of starvation or high LDL exposure, the lipid droplet core can undergo an amorphous to liquid-crystalline phase transformation. Herein, we study the consequences of this transformation on the physical properties of lipid droplets that are thought to regulate protein association. Using simulations of different sterol-ester concentrations, we have captured the liquid-crystalline phase transformation at the molecular level, highlighting the alignment of sterol esters in alternating orientations to form concentric layers. We demonstrate how ordering in the core permeates into the neutral lipid/phospholipid interface, changing the magnitude and nature of neutral lipid intercalation and inducing ordering in the phospholipid monolayer. Increased phospholipid packing is concomitant with altered surface properties, including smaller area per phospholipid and substantially reduced packing defects. Additionally, the ordering of sterol esters in the core causes less hydration in more ordered regions. We discuss these findings in the context of their expected consequences for preferential protein recruitment to lipid droplets under different metabolic conditions.

Computational Studies of Lipid Droplets

Lipid droplets (LDs) are intracellular organelles whose primary function is energy storage. Known to emerge from the endoplasmic reticulum (ER) bilayer, LDs have a unique structure with a core consisting of neutral lipids, triacylglycerol (TG) or sterol esters (SE), surrounded by a phospholipid (PL) monolayer and decorated by proteins that come and go throughout their complex lifecycle. In this Feature Article, we review recent developments in computational studies of LDs, a rapidly growing area of research. We highlight how molecular dynamics (MD) simulations have provided valuable molecular-level insight into LD targeting and LD biogenesis. Additionally, we review the physical properties of TG from different force fields compared with experimental data. Possible future directions and challenges are discussed.

The CYTOLD and ERTOLD pathways for lipid droplet-protein targeting

Lipid droplets (LDs) are the main organelles for lipid storage, and their surfaces contain unique proteins with diverse functions, including those that facilitate the deposition and mobilization of LD lipids. Among organelles, LDs have an unusual structure with an organic, hydrophobic oil phase covered by a phospholipid monolayer. The unique properties of LD monolayer surfaces require proteins to localize to LDs by distinct mechanisms. Here we review the two pathways known to mediate direct LD protein localization: the CYTOLD pathway mediates protein targeting from the cytosol toLDs, and the ERTOLD pathway functions in protein targeting from the endoplasmic reticulum toLDs. We describe the emerging principles for each targeting pathway in animal cells and highlight open questions in the field.

Toward a Multipathway Perspective: pH-Dependent Kinetic Selection of Competing Pathways and the Role of the Internal Glutamate in Cl-/H+ Antiporters

Canonical descriptions of multistep biomolecular transformations generally follow a single-pathway viewpoint, with a series of transitions through intermediates converting reactants to products or repeating a conformational cycle. However, mounting evidence suggests that more complexity and pathway heterogeneity are mechanistically relevant due to the statistical distribution of multiple interconnected rate processes. Making sense of such pathway complexity remains a significant challenge. To better understand the role and relevance of pathway heterogeneity, we herein probe the chemical reaction network of a Cl-/H+ antiporter, ClC-ec1, and analyze reaction pathways using multiscale kinetic modeling (MKM). This approach allows us to describe the nature of the competing pathways and how they change as a function of pH. We reveal that although pH-dependent Cl-/H+ transport rates are largely regulated by the charge state of amino acid E148, the charge state of E203 determines relative contributions from coexisting pathways and can shift the flux pH-dependence. The selection of pathways via E203 explains how ionizable mutations (D/H/K/R) would impact the ClC-ec1 bioactivity from a kinetic perspective and lends further support to the indispensability of an internal glutamate in ClC antiporters. Our results demonstrate how quantifying the kinetic selection of competing pathways under varying conditions leads to a deeper understanding of the Cl-/H+ exchange mechanism and can suggest new approaches for mechanistic control.

Improving the accuracy and convergence of drug permeation simulations via machine-learned collective variables

Understanding the permeation of biomolecules through cellular membranes is critical for many biotechnological applications, including targeted drug delivery, pathogen detection, and the development of new antibiotics. To this end, computer simulations are routinely used to probe the underlying mechanisms of membrane permeation. Despite great progress and continued development, permeation simulations of realistic systems (e.g., more complex drug molecules or biologics through heterogeneous membranes) remain extremely challenging if not intractable. In this work, we combine molecular dynamics simulations with transition-tempered metadynamics and techniques from the variational approach to conformational dynamics to study the permeation mechanism of a drug molecule, trimethoprim, through a multicomponent membrane. We show that collective variables (CVs) obtained from an unsupervised machine learning algorithm called time-structure based Independent Component Analysis (tICA) improve performance and substantially accelerate convergence of permeation potential of mean force (PMF) calculations. The addition of cholesterol to the lipid bilayer is shown to increase both the width and height of the free energy barrier due to a condensing effect (lower area per lipid) and increase bilayer thickness. Additionally, the tICA CVs reveal a subtle effect of cholesterol increasing the resistance to permeation in the lipid head group region, which is not observed when canonical CVs are used. We conclude that the use of tICA CVs can enable more efficient PMF calculations with additional insight into the permeation mechanism.

Stressed Lipid Droplets: How Neutral Lipids Relieve Surface Tension and Membrane Expansion Drives Protein Association

Lipid droplets (LDs) are intracellular storage organelles composed of neutral lipids, such as triacylglycerol (TG), surrounded by a phospholipid (PL) monolayer decorated with specific proteins. Herein, we investigate the mechanism of protein association during LD and bilayer membrane expansion. We find that the neutral lipids play a dynamic role in LD expansion by further intercalating with the PL monolayer to create more surface-oriented TG molecules (SURF-TG). This interplay both reduces high surface tension incurred during LD budding or growth and also creates expansion-specific surface features for protein recognition. We then show that the autoinhibitory (AI) helix of CTP:phosphocholine cytidylyltransferase, a protein known to target expanding monolayers and bilayers, preferentially associates with large packing defects in a sequence-specific manner. Despite the presence of three phenylalanines, the initial binding with bilayers is predominantly mediated by the sole tryptophan due to its preference for membrane interfaces. Subsequent association is dependent on the availability of large, neighboring defects that can accommodate the phenylalanines, which are more probable in the stressed systems. Tryptophan, once fully associated, preferentially interacts with the glycerol moiety of SURF-TG in LDs. The calculation of AI binding free energy, hydrogen bonding and depth analysis, and in silico mutation experiments support the findings. Hence, SURF-TG can both reduce surface tension and mediate protein association, facilitating class II protein recruitment during LD expansion.

The Surface and Hydration Properties of Lipid Droplets

Lipid droplets (LDs) are energy storage organelles composed of neutral lipids, such as triacylglycerol (TG) and sterol esters, surrounded by a phospholipid (PL) monolayer. Their central role in metabolism, complex life cycle, and unique lipid monolayer surface have garnered great attention over the last decade. In this article, results from the largest and longest all-atom simulations to date suggest that 5-8% of the LD surface is occupied by TG molecules, a number that exceeds the maximal solubility reported for TGs in PL bilayers (2.8%). Two distinct classes of TG molecules that interact with the LD monolayer are found. Those at the monolayer surface (SURF-TG) are ordered like PLs with the glycerol moiety exposed to water, creating a significant amount of chemically unique packing defects, and the acyl chains extended toward the LD center. In contrast, the TGs that intercalate just into the PL tail region (CORE-TG) are disordered and increase the amount of PL packing defects and the PL tail order. The degree of interdigitation caused by CORE-TG is stable and determines the width of the TG-PL overlap, whereas that caused by SURF-TG fluctuates and is highly correlated with the area per PL or the expansion of the monolayer. Thus, when the supply of PLs is limited, SURF-TG may reduce surface tension by behaving as a secondary membrane component. The hydration properties of the simulated LD systems demonstrate ∼10 times more water in the LD core than previously reported. Collectively, the reported surface and hydration properties are expected to play a direct role in the mechanisms by which proteins target and interact with LDs.

Understanding and Tracking the Excess Proton in Ab Initio Simulations; Insights from IR Spectra

Proton transport in aqueous media is ubiquitously important in chemical and biological processes. Although ab initio molecular dynamics (AIMD) simulations have made great progress in characterizing proton transport, there has been a long-standing challenge in defining and tracking the excess proton, or more properly, the center of excess charge (CEC) created when a hydrogen nucleus distorts the electron distributions of water molecules in a delocalized and highly dynamic nature. Yet, defining (and biasing) such a CEC is essential when combining AIMD with enhanced sampling methods to calculate the relevant macroscopic properties via free-energy landscapes, which is the standard practice for most processes of interest. Several CEC formulas have been proposed and used, but none have yet been systematically tested or rigorously derived. In this paper, we show that the CEC can be used as a computational tool to disentangle IR features of the solvated excess proton from its surrounding solvent, and in turn, how correlating the features in the excess charge spectrum with the behavior of CEC in simulations enables a systematic evaluation of various CEC definitions. We present a new definition of CEC and show how it overcomes the limitations of those currently available both from a spectroscopic point of view and from a practical perspective of performance in enhanced sampling simulations.

Local conformational dynamics regulating transport properties of a Cl- /H+ antiporter

ClC-ec1 is a Cl- /H+ antiporter that exchanges Cl- and H+ ions across the membrane. Experiments have demonstrated that several mutations, including I109F, decrease the Cl- and H+ transport rates by an order of magnitude. Using reactive molecular dynamics simulations of explicit proton transport across the central region in the I109F mutant, a two-dimensional free energy profile has been constructed that is consistent with the experimental transport rates. The importance of a phenylalanine gate formed by F109 and F357 and its influence on hydration connectivity through the central proton transport pathway is revealed. This work demonstrates how seemingly subtle changes in local conformational dynamics can dictate hydration changes and thus transport properties. © 2019 Wiley Periodicals, Inc.

Dynamic Protonation Dramatically Affects the Membrane Permeability of Drug-like Molecules

Permeability (P_m) across biological membranes is of fundamental importance and a key factor in drug absorption, distribution, and development. Although the majority of drugs will be charged at some point during oral delivery, our understanding of membrane permeation by charged species is limited. The canonical model assumes that only neutral molecules partition into and passively permeate across membranes, but there is mounting evidence that these processes are also facile for certain charged species. However, it is unknown whether such ionizable permeants dynamically neutralize at the membrane surface or permeate in their charged form. To probe protonation-coupled permeation in atomic detail, we herein apply continuous constant-pH molecular dynamics along with free energy sampling to study the permeation of a weak base propranolol (PPL), and evaluate the impact of including dynamic protonation on P_m. The simulations reveal that PPL dynamically neutralizes at the lipid-tail interface, which dramatically influences the permeation free energy landscape and explains why the conventional model overestimates the assigned intrinsic permeability. We demonstrate how fixed-charge-state simulations can account for this effect, and propose a revised model that better describes pH-coupled partitioning and permeation. Our results demonstrate how dynamic changes in protonation state may play a critical role in the permeation of ionizable molecules, including pharmaceuticals and drug-like molecules, thus requiring a revision of the standard picture.

Proton-Induced Conformational and Hydration Dynamics in the Influenza A M2 Channel

The influenza A M2 protein is an acid-activated proton channel responsible for acidification of the inside of the virus, a critical step in the viral life cycle. This channel has four central histidine residues that form an acid-activated gate, binding protons from the outside until an activated state allows proton transport to the inside. While previous work has focused on proton transport through the channel, the structural and dynamic changes that accompany proton flux and enable activation have yet to be resolved. In this study, extensive Multiscale Reactive Molecular Dynamics simulations with explicit Grotthuss-shuttling hydrated excess protons are used to explore detailed molecular-level interactions that accompany proton transport in the +0, + 1, and +2 histidine charge states. The results demonstrate how the hydrated excess proton strongly influences both the protein and water hydrogen-bonding network throughout the channel, providing further insight into the channel's acid-activation mechanism and rectification behavior. We find that the excess proton dynamically, as a function of location, shifts the protein structure away from its equilibrium distributions uniquely for different pH conditions consistent with acid-activation. The proton distribution in the xy-plane is also shown to be asymmetric about the channel's main axis, which has potentially important implications for the mechanism of proton conduction and future drug design efforts.

Mycolactone Toxin Membrane Permeation: Atomistic versus Coarse-Grained MARTINI Simulations

Mycolactone, a cytotoxic and immunosuppressive macrolide produced by Mycobacterium ulcerans, is the central virulent factor in the skin disease Buruli ulcer. This multifunctional cytotoxin affects fundamental cellular processes such as cell adhesion, immune response, and cell death by targeting various cellular structures. Developing effective diagnostics that target mycolactone has been challenging, potentially because of suspected interactions with lipophilic architectures, including membranes. To better understand the pathogenesis of Buruli ulcer disease, aid in the development of diagnostics, and learn how amphiphiles in general use lipid trafficking to navigate the host environment, we seek to understand the nature of mycolactone-membrane interactions. Herein, we characterize how the two dominant isomers of mycolactone (A and B) interact with and permeate DPPC membranes with all-atom molecular dynamics simulations employing transition-tempered metadynamics and compare these results to those obtained by MARTINI coarse-grained simulations. Our all-atom simulations reveal that both isomers have a strong preference to associate with the membrane, although their mechanisms and energetics of membrane permeation differ slightly. Water molecules are found to play an important role in the permeation process. Although the MARTINI coarse-grained simulations give the correct free energy of membrane association, they fail to capture the mechanism of permeation and role of water during permeation as seen in all-atom simulations.

Modulating the Chemical Transport Properties of a Transmembrane Antiporter via Alternative Anion Flux

ClC-ec1 is a prototype of the ClC antiporters, proteins that stoichiometrically exchange Cl^– and H⁺ ions in opposite directions across a membrane. It has been shown that other polyatomic anions, such as NO₃^– and SCN^–, can also be transported by ClC-ec1, but with partially or completely uncoupled proton flux. Herein, with the help of multiscale computer simulations in which the Grotthuss mechanism of proton transport (PT) is treated explicitly, we demonstrate how the chemical nature of these anions alters the coupling mechanism and qualitatively explain the shifts in the ion stoichiometry. Multidimensional free energy profiles for PT and the coupled changes in hydration are presented for NO₃^– and SCN^–. The calculated proton conductances agree with experiment, showing reduced or abolished proton flux. Surprisingly, the proton affinity of the anion is less influential on the PT, while its size and interactions with the protein significantly alter hydration and shift its influence on PT from facilitating to inhibiting. We find that the hydration of the cavity below the anion is relatively fast, but connecting the water network past the steric hindrance of these polyatomic anions is quite slow. Hence, the most relevant coordinate to the PT free energy barrier is the water connectivity along the PT pathway, but importantly only in the presence of the excess proton, and this coordinate is significantly affected by the nature of the bound anion. This work again demonstrates how PT is intrinsically coupled with protein cavity hydration changes as well as influenced by the protein environment. It additionally suggests ways in which ion exchange can be modulated and exchange stoichiometries altered.

Molecular transport through membranes: Accurate permeability coefficients from multidimensional potentials of mean force and local diffusion constants

Estimating the permeability coefficient of small molecules through lipid bilayer membranes plays an important role in the development of effective drug candidates. In silico simulations can produce acceptable relative permeability coefficients for a series of small molecules; however, the absolute permeability coefficients from simulations are usually off by orders of magnitude. In addition to differences between the lipid bilayers used in vitro and in silico, the poor convergence of permeation free energy profiles and over-simplified diffusion models have contributed to these discrepancies. In this paper, we present a multidimensional inhomogeneous solubility-diffusion model to study the permeability of a small molecule drug (trimethoprim) passing through a POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) lipid bilayer. Our approach improves the permeation model in three ways: First, the free energy profile (potential of mean force, PMF) is two-dimensional in two key coordinates rather than simply one-dimensional along the direction normal to the bilayer. Second, the 2-D PMF calculation has improved convergence due to application of the recently developed transition-tempered metadynamics with randomly initialized replicas, while third, the local diffusivity coefficient was calculated along the direction of the minimum free energy path on the two-dimensional PMF. The permeability is then calculated as a line integral along the minimum free energy path of the PMF. With this approach, we report a considerably more accurate permeability coefficient (only 2-5 times larger than the experimental value). We also compare our approach with the common practice of computing permeability coefficients based only on the translation of the center of mass of the drug molecule. Our paper concludes with a discussion of approaches for minimizing the computational cost for the purpose of more rapidly screening a large number of drug candidate molecules.

Multiscale Kinetic Modeling Reveals an Ensemble of Cl-/H+ Exchange Pathways in ClC-ec1 Antiporter

Despite several years of research, the ion exchange mechanisms in chloride/proton antiporters and many other coupled transporters are not yet understood at the molecular level. Here, we present a novel approach to kinetic modeling and apply it to ion exchange in ClC-ec1. Our multiscale kinetic model is developed by (1) calculating the state-to-state rate coefficients with reactive and polarizable molecular dynamics simulations, (2) optimizing these rates in a global kinetic network, and (3) predicting new electrophysiological results. The model shows that the robust Cl:H exchange ratio (2.2:1) can indeed arise from kinetic coupling without large protein conformational changes, indicating a possible facile evolutionary connection to chloride channels. The E148 amino acid residue is shown to couple chloride and proton transport through protonation-dependent blockage of the central anion binding site and an anion-dependent pK_a value, which influences proton transport. The results demonstrate how an ensemble of different exchange pathways, as opposed to a single series of transitions, culminates in the macroscopic observables of the antiporter, such as transport rates, chloride/proton stoichiometry, and pH dependence.

Understanding the essential proton-pumping kinetic gates and decoupling mutations in cytochrome c oxidase

Cytochrome c oxidase (CcO) catalyzes the reduction of oxygen to water and uses the released free energy to pump protons against the transmembrane proton gradient. To better understand the proton-pumping mechanism of the wild-type (WT) CcO, much attention has been given to the mutation of amino acid residues along the proton translocating D-channel that impair, and sometimes decouple, proton pumping from the chemical catalysis. Although their influence has been clearly demonstrated experimentally, the underlying molecular mechanisms of these mutants remain unknown. In this work, we report multiscale reactive molecular dynamics simulations that characterize the free-energy profiles of explicit proton transport through several important D-channel mutants. Our results elucidate the mechanisms by which proton pumping is impaired, thus revealing key kinetic gating features in CcO. In the N139T and N139C mutants, proton back leakage through the D-channel is kinetically favored over proton pumping due to the loss of a kinetic gate in the N139 region. In the N139L mutant, the bulky L139 side chain inhibits timely reprotonation of E286 through the D-channel, which impairs both proton pumping and the chemical reaction. In the S200V/S201V double mutant, the proton affinity of E286 is increased, which slows down both proton pumping and the chemical catalysis. This work thus not only provides insight into the decoupling mechanisms of CcO mutants, but also explains how kinetic gating in the D-channel is imperative to achieving high proton-pumping efficiency in the WT CcO.

Acid activation mechanism of the influenza A M2 proton channel

The homotetrameric influenza A M2 channel (AM2) is an acid-activated proton channel responsible for the acidification of the influenza virus interior, an important step in the viral lifecycle. Four histidine residues (His37) in the center of the channel act as a pH sensor and proton selectivity filter. Despite intense study, the pH-dependent activation mechanism of the AM2 channel has to date not been completely understood at a molecular level. Herein we have used multiscale computer simulations to characterize (with explicit proton transport free energy profiles and their associated calculated conductances) the activation mechanism of AM2. All proton transfer steps involved in proton diffusion through the channel, including the protonation/deprotonation of His37, are explicitly considered using classical, quantum, and reactive molecular dynamics methods. The asymmetry of the proton transport free energy profile under high-pH conditions qualitatively explains the rectification behavior of AM2 (i.e., why the inward proton flux is allowed when the pH is low in viral exterior and high in viral interior, but outward proton flux is prohibited when the pH gradient is reversed). Also, in agreement with electrophysiological results, our simulations indicate that the C-terminal amphipathic helix does not significantly change the proton conduction mechanism in the AM2 transmembrane domain; the four transmembrane helices flanking the channel lumen alone seem to determine the proton conduction mechanism.

The Origin of Coupled Chloride and Proton Transport in a Cl-/H+ Antiporter

The ClC family of transmembrane proteins functions throughout nature to control the transport of Cl^– ions across biological membranes. ClC-ec1 from Escherichia coli is an antiporter, coupling the transport of Cl^– and H⁺ ions in opposite directions and driven by the concentration gradients of the ions. Despite keen interest in this protein, the molecular mechanism of the Cl^–/H⁺ coupling has not been fully elucidated. Here, we have used multiscale simulation to help identify the essential mechanism of the Cl^–/H⁺ coupling. We find that the highest barrier for proton transport (PT) from the intra- to extracellular solution is attributable to a chemical reaction, the deprotonation of glutamic acid 148 (E148). This barrier is significantly reduced by the binding of Cl^– in the “central” site (Cl^–_cen), which displaces E148 and thereby facilitates its deprotonation. Conversely, in the absence of Cl^–_cen E148 favors the “down” conformation, which results in a much higher cumulative rotation and deprotonation barrier that effectively blocks PT to the extracellular solution. Thus, the rotation of E148 plays a critical role in defining the Cl^–/H⁺ coupling. As a control, we have also simulated PT in the ClC-ec1 E148A mutant to further understand the role of this residue. Replacement with a non-protonatable residue greatly increases the free energy barrier for PT from E203 to the extracellular solution, explaining the experimental result that PT in E148A is blocked whether or not Cl^–_cen is present. The results presented here suggest both how a chemical reaction can control the rate of PT and also how it can provide a mechanism for a coupling of the two ion transport processes.

Computationally Efficient Multiscale Reactive Molecular Dynamics to Describe Amino Acid Deprotonation in Proteins

An important challenge in the simulation of biomolecular systems is a quantitative description of the protonation and deprotonation process of amino acid residues. Despite the seeming simplicity of adding or removing a positively charged hydrogen nucleus, simulating the actual protonation/deprotonation process is inherently difficult. It requires both the explicit treatment of the excess proton, including its charge defect delocalization and Grotthuss shuttling through inhomogeneous moieties (water and amino residues), and extensive sampling of coupled condensed phase motions. In a recent paper (J. Chem. Theory Comput. 2014, 10, 2729-2737), a multiscale approach was developed to map high-level quantum mechanics/molecular mechanics (QM/MM) data into a multiscale reactive molecular dynamics (MS-RMD) model in order to describe amino acid deprotonation in bulk water. In this article, we extend the fitting approach (called FitRMD) to create MS-RMD models for ionizable amino acids within proteins. The resulting models are shown to faithfully reproduce the free energy profiles of the reference QM/MM Hamiltonian for PT inside an example protein, the ClC-ec1 H(+)/Cl(-) antiporter. Moreover, we show that the resulting MS-RMD models are computationally efficient enough to then characterize more complex 2-dimensional free energy surfaces due to slow degrees of freedom such as water hydration of internal protein cavities that can be inherently coupled to the excess proton charge translocation. The FitRMD method is thus shown to be an effective way to map ab initio level accuracy into a much more computationally efficient reactive MD method in order to explicitly simulate and quantitatively describe amino acid protonation/deprotonation in proteins.

Optimized Radii for Poisson-Boltzmann Calculations with the AMBER Force Field

Implicit solvent models are a standard tool for assessing the electrostatics of biomolecular systems. The accuracy of quantitative predictions, such as pKa values, transfer free energies, binding energies, and solvation forces, is strongly dependent on one's choice of continuum parameters:  the solute charges, dielectric coefficient, and radii, which define the dielectric boundary. To ensure quantitative accuracy, these parameters can be benchmarked against explicit solvent simulations. Here we present two sets of optimized radii to define either abrupt or cubic-spline smoothed dielectric boundaries in Poisson-Boltzmann calculations of protein systems with AMBER (parm99) charges. Spline smoothing stabilizes the electrostatic potential at the molecular surface, allowing for continuum force calculations. Most implementations, however, require significantly different radii than the abrupt boundary surfaces. The optimal continuum radii are initially approximated from the solvent radial charge distribution surrounding each atom type. A genetic algorithm is then used to fine-tune the starting values to reproduce charging free energies measured from explicit solvent simulations. The optimized radii are tested on four protein-like polypeptides. The results show increased accuracy of molecular solvation energies and atomic forces relative to commonly used continuum parameter sets. These radii are suitable for Poisson-Boltzmann calculations with the AMBER force field and offer energetic congruence to any model that combines molecular mechanics and Poisson-Boltzmann solvation energies.

Optimizing the Poisson Dielectric Boundary with Explicit Solvent Forces and Energies: Lessons Learned with Atom-Centered Dielectric Functions

Accurate implicit solvent models require parameters that have been optimized in a system- or atom-specific manner on the basis of experimental data or more rigorous explicit solvent simulations. Models based on the Poisson or Poisson-Boltzmann equation are particularly sensitive to the nature and location of the boundary which separates the low dielectric solute from the high dielectric solvent. Here, we present a novel method for optimizing the solute radii, which define the dielectric boundary, on the basis of forces and energies from explicit solvent simulations. We use this method to optimize radii for protein systems defined by AMBER ff99 partial charges and a spline-smoothed solute surface. The spline-smoothed surface is an atom-centered dielectric function that enables stable and efficient force calculations. We explore the relative performance of radii optimized with forces alone and those optimized with forces and energies. We show that our radii reproduce the explicit solvent forces and energies more accurately than four other parameter sets commonly used in conjunction with the AMBER force field, each of which has been appropriately scaled for spline-smoothed surfaces. Finally, we demonstrate that spline-smoothed surfaces show surprising accuracy for small, compact systems but may have limitations for highly solvated protein systems. The optimization method presented here is efficient and applicable to any system with explicit solvent parameters. It can be used to determine the optimal continuum parameters when experimental solvation energies are unavailable and the computational costs of explicit solvent charging free energies are prohibitive.

Hydrated Excess Protons Can Create Their Own Water Wires

Grotthuss shuttling of an excess proton charge defect through hydrogen bonded water networks has long been the focus of theoretical and experimental studies. In this work we show that there is a related process in which water molecules move ("shuttle") through a hydrated excess proton charge defect in order to wet the path ahead for subsequent proton charge migration. This process is illustrated through reactive molecular dynamics simulations of proton transport through a hydrophobic nanotube, which penetrates through a hydrophobic region. Surprisingly, before the proton enters the nanotube, it starts "shooting" water molecules into the otherwise dry space via Grotthuss shuttling, effectively creating its own water wire where none existed before. As the proton enters the nanotube (by 2-3 Å), it completes the solvation process, transitioning the nanotube to the fully wet state. By contrast, other monatomic cations (e.g., K(+)) have just the opposite effect, by blocking the wetting process and making the nanotube even drier. As the dry nanotube gradually becomes wet when the proton charge defect enters it, the free energy barrier of proton permeation through the tube via Grotthuss shuttling drops significantly. This finding suggests that an important wetting mechanism may influence proton translocation in biological systems, i.e., one in which protons "create" their own water structures (water "wires") in hydrophobic spaces (e.g., protein pores) before migrating through them. An existing water wire, e.g., one seen in an X-ray crystal structure or MD simulations without an explicit excess proton, is therefore not a requirement for protons to transport through hydrophobic spaces.

Multiscale Reactive Molecular Dynamics for Absolute pKa Predictions and Amino Acid Deprotonation

Accurately calculating a weak acid's pK_a from simulations remains a challenging task. We report a multiscale theoretical approach to calculate the free energy profile for acid ionization, resulting in accurate absolute pK_a values in addition to insights into the underlying mechanism. Importantly, our approach minimizes empiricism by mapping electronic structure data (QM/MM forces) into a reactive molecular dynamics model capable of extensive sampling. Consequently, the bulk property of interest (the absolute pK_a) is the natural consequence of the model, not a parameter used to fit it. This approach is applied to create reactive models of aspartic and glutamic acids. We show that these models predict the correct pK_a values and provide ample statistics to probe the molecular mechanism of dissociation. This analysis shows changes in the solvation structure and Zundel-dominated transitions between the protonated acid, contact ion pair, and bulk solvated excess proton.

Benchmark Study of the SCC-DFTB Approach for a Biomolecular Proton Channel

The self-consistent charge density functional tight binding (SCC-DFTB) method has been increasingly applied to study proton transport (PT) in biological environments. However, recent studies revealing some significant limitations of SCC-DFTB for proton and hydroxide solvation and transport in bulk aqueous systems call into question its accuracy for simulating PT in biological systems. The current work benchmarks the SCC-DFTB/MM method against more accurate DFT/MM by simulating PT in a synthetic leucine-serine channel (LS2), which emulates the structure and function of biomolecular proton channels. It is observed that SCC-DFTB/MM produces over-coordinated and less structured pore water, an over-coordinated excess proton, weak hydrogen bonds around the excess proton charge defect and qualitatively different PT dynamics. Similar issues are demonstrated for PT in a carbon nanotube, indicating that the inaccuracies found for SCC-DFTB are not due to the point charge based QM/MM electrostatic coupling scheme, but rather to the approximations of the semiempirical method itself. The results presented in this work highlight the limitations of the present form of the SCC-DFTB/MM approach for simulating PT processes in biological protein or channel-like environments, while providing benchmark results that may lead to an improvement of the underlying method.

Intricate role of water in proton transport through cytochrome c oxidase

Cytochrome c oxidase (CytcO), the final electron acceptor in the respiratory chain, catalyzes the reduction of O(2) to H(2)O while simultaneously pumping protons across the inner mitochondrial or bacterial membrane to maintain a transmembrane electrochemical gradient that drives, for example, ATP synthesis. In this work mutations that were predicted to alter proton translocation and enzyme activity in preliminary computational studies are characterized with extensive experimental and computational analysis. The mutations were introduced in the D pathway, one of two proton-uptake pathways, in CytcO from Rhodobacter sphaeroides . Serine residues 200 and 201, which are hydrogen-bonded to crystallographically resolved water molecules halfway up the D pathway, were replaced by more bulky hydrophobic residues (Ser200Ile, Ser200Val/Ser201Val, and Ser200Val/Ser201Tyr) to query the effects of changing the local structure on enzyme activity as well as proton uptake, release, and intermediate transitions. In addition, the effects of these mutations on internal proton transfer were investigated by blocking proton uptake at the pathway entrance (Asp132Asn replacement in addition to the above-mentioned mutations). Even though the overall activities of all mutant CytcO's were lowered, both the Ser200Ile and Ser200Val/Ser201Val variants maintained the ability to pump protons. The lowered activities were shown to be due to slowed oxidation kinetics during the P(R) → F and F → O transitions (P(R) is the "peroxy" intermediate formed at the catalytic site upon reaction of the four-electron-reduced CytcO with O(2), F is the oxoferryl intermediate, and O is the fully oxidized CytcO). Furthermore, the P(R) → F transition is shown to be essentially pH independent up to pH 12 (i.e., the apparent pK(a) of Glu286 is increased from 9.4 by at least 3 pK(a) units) in the Ser200Val/Ser201Val mutant. Explicit simulations of proton transport in the mutated enzymes revealed that the solvation dynamics can cause intriguing energetic consequences and hence provide mechanistic insights that would never be detected in static structures or simulations of the system with fixed protonation states (i.e., lacking explicit proton transport). The results are discussed in terms of the proton-pumping mechanism of CytcO.

Role of charge transfer in the structure and dynamics of the hydrated proton

Although it has long been recognized that multiple water molecules strongly associate with an extra proton in bulk water, some models and conceptual frameworks continue to utilize the classical hydronium ion (H(3)O(+)) as a fundamental building block. In this work, the nature of the hydronium ion in aqueous systems is examined using an ab initio energy decomposition analysis (EDA) that evaluates both the magnitude of and energetic stabilization due to charge transfer among H(3)O(+) and the surrounding water molecules. The EDA is performed on structures extracted from dynamical bulk-phase simulations and used to determine how frequently the pure hydronium ion, where the excess charge is primarily localized on H(3)O(+), occurs under dynamic conditions. The answer is essentially never. The energetic stabilization of H(3)O(+) due to charge delocalization to neighboring water molecules is found to be much larger (16-49 kcal/mol) than for other ions (even Li(+)) and to constitute a substantial portion (20-52%) of the complex's binding energy. The charge defect is also shown to have intrinsic dynamical asymmetry and to display dynamical signatures that can be related to features appearing in IR spectra.

Coarse-graining in interaction space: a systematic approach for replacing long-range electrostatics with short-range potentials

A short-range effective potential for long-range electrostatic interactions in homogeneously disordered condensed phase systems has been determined with a novel approach to coarse-graining in interaction space. As opposed to coarse-graining the system resolution, this approach "coarsens" the system's interactions by mapping multiple configurations of an accurate long-range atomistic potential onto a more efficient, short-range effective potential with a force-matching (FM) method. Developing an empirical potential in this manner is fundamentally different from existing strategies because it utilizes condensed-phase (as opposed to gas-phase) atomistic interactions to determine general pair potentials defined on distance meshes (as opposed to fitting predetermined functional forms). The resulting short-range ( approximately 10 A) effective potential reproduces structural, dynamical, and many thermodynamic properties of liquid water, ions in water, and hydrophobes in water, with unprecedented accuracy. The effective potential is also shown to be transferable to a nonaqueous molten salt system. With continued development, such effective potentials may provide an accurate and highly efficient alternative to Ewald-based long-range electrostatics methods.

Limitations of atom-centered dielectric functions in implicit solvent models

Many recent advances in Poisson-Boltzmann and generalized Born implicit solvent models have used atom-centered polynomial or Gaussian functions to define the boundary separating low and high dielectric regions. In contrast to the Lee and Richards molecular surface, atom-centered surfaces result in interatomic crevices and buried pockets of high dielectric which are too small for a solvent molecule to occupy. We show that these interstitial high dielectric regions are of significant magnitude in globular proteins, that they artificially increase solvation energies, and that they distort the free energy surface of nonbonded interactions. These results suggest that implicit solvent dielectric functions must exclude interstitial high dielectric regions in order to yield physically meaningful results.

Proton solvation and transport in aqueous and biomolecular systems: insights from computer simulations

The excess proton in aqueous media plays a pivotal role in many fundamental chemical (e.g., acid-base chemistry) and biological (e.g., bioenergetics and enzyme catalysis) processes. Understanding the hydrated proton is, therefore, crucial for chemistry, biology, and materials sciences. Although well studied for over 200 years, excess proton solvation and transport remains to this day mysterious, surprising, and perhaps even misunderstood. In this feature article, various efforts to address this problem through computer modeling and simulation will be described. Applications of computer simulations to a number of important and interesting systems will be presented, highlighting the roles of charge delocalization and Grotthuss shuttling, a phenomenon unique in many ways to the excess proton in water.

Coupling hydrophobicity, dispersion, and electrostatics in continuum solvent models

An implicit solvent model is presented that couples hydrophobic, dispersion, and electrostatic solvation energies by minimizing the system Gibbs free energy with respect to the solvent volume exclusion function. The solvent accessible surface is the output of the theory. The method is illustrated with the solvation of simple solutes on different length scales and captures the sensitivity of hydration to the particular form of the solute-solvent interactions in agreement with recent computer simulations.

Coupling nonpolar and polar solvation free energies in implicit solvent models

Recent studies on the solvation of atomistic and nanoscale solutes indicate that a strong coupling exists between the hydrophobic, dispersion, and electrostatic contributions to the solvation free energy, a facet not considered in current implicit solvent models. We suggest a theoretical formalism which accounts for coupling by minimizing the Gibbs free energy of the solvent with respect to a solvent volume exclusion function. The resulting differential equation is similar to the Laplace-Young equation for the geometrical description of capillary interfaces but is extended to microscopic scales by explicitly considering curvature corrections as well as dispersion and electrostatic contributions. Unlike existing implicit solvent approaches, the solvent accessible surface is an output of our model. The presented formalism is illustrated on spherically or cylindrically symmetrical systems of neutral or charged solutes on different length scales. The results are in agreement with computer simulations and, most importantly, demonstrate that our method captures the strong sensitivity of solvent expulsion and dewetting to the particular form of the solvent-solute interactions.

Revisiting free energy calculations: a theoretical connection to MM/PBSA and direct calculation of the association free energy

The prediction of absolute ligand-receptor binding affinities is essential in a wide range of biophysical queries, from the study of protein-protein interactions to structure-based drug design. End-point free energy methods, such as the Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) model, have received much attention and widespread application in recent literature. These methods benefit from computational efficiency as only the initial and final states of the system are evaluated, yet there remains a need for strengthening their theoretical foundation. Here a clear connection between statistical thermodynamics and end-point free energy models is presented. The importance of the association free energy, arising from one molecule's loss of translational and rotational freedom from the standard state concentration, is addressed. A novel method for calculating this quantity directly from a molecular dynamics simulation is described. The challenges of accounting for changes in the protein conformation and its fluctuations from separate simulations are discussed. A simple first-order approximation of the configuration integral is presented to lay the groundwork for future efforts. This model has been applied to FKBP12, a small immunophilin that has been widely studied in the drug industry for its potential immunosuppressive and neuroregenerative effects.