Clustering Protein Binding Pockets and Identifying Potential Drug Interactions: A Novel Ligand-Based Featurization Method.

Protein-ligand interactions are essential to drug discovery and drug development efforts. Desirable on-target or multitarget interactions are the first step in finding an effective therapeutic, while undesirable off-target interactions are the first step in assessing safety. In this work, we introduce a novel ligand-based featurization and mapping of human protein pockets to identify closely related protein targets and to project novel drugs into a hybrid protein-ligand feature space to identify their likely protein interactions. Using structure-based template matches from PDB, protein pockets are featured by the ligands that bind to their best co-complex template matches. The simplicity and interpretability of this approach provide a granular characterization of the human proteome at the protein-pocket level instead of the traditional protein-level characterization by family, function, or pathway. We demonstrate the power of this featurization method by clustering a subset of the human proteome and evaluating the predicted cluster associations of over 7000 compounds.

Improved Protein-Ligand Binding Affinity Prediction with Structure-Based Deep Fusion Inference.

Predicting accurate protein-ligand binding affinities is an important task in drug discovery but remains a challenge even with computationally expensive biophysics-based energy scoring methods and state-of-the-art deep learning approaches. Despite the recent advances in the application of deep convolutional and graph neural network-based approaches, it remains unclear what the relative advantages of each approach are and how they compare with physics-based methodologies that have found more mainstream success in virtual screening pipelines. We present fusion models that combine features and inference from complementary representations to improve binding affinity prediction. This, to our knowledge, is the first comprehensive study that uses a common series of evaluations to directly compare the performance of three-dimensional (3D)-convolutional neural networks (3D-CNNs), spatial graph neural networks (SG-CNNs), and their fusion. We use temporal and structure-based splits to assess performance on novel protein targets. To test the practical applicability of our models, we examine their performance in cases that assume that the crystal structure is not available. In these cases, binding free energies are predicted using docking pose coordinates as the inputs to each model. In addition, we compare these deep learning approaches to predictions based on docking scores and molecular mechanic/generalized Born surface area (MM/GBSA) calculations. Our results show that the fusion models make more accurate predictions than their constituent neural network models as well as docking scoring and MM/GBSA rescoring, with the benefit of greater computational efficiency than the MM/GBSA method. Finally, we provide the code to reproduce our results and the parameter files of the trained models used in this work. The software is available as open source at https://github.com/llnl/fast. Model parameter files are available at ftp://gdo-bioinformatics.ucllnl.org/fast/pdbbind2016_model_checkpoints/.

Predicting Small Molecule Transfer Free Energies by Combining Molecular Dynamics Simulations and Deep Learning.

Accurately predicting small molecule partitioning and hydrophobicity is critical in the drug discovery process. There are many heterogeneous chemical environments within a cell and entire human body. For example, drugs must be able to cross the hydrophobic cellular membrane to reach their intracellular targets, and hydrophobicity is an important driving force for drug-protein binding. Atomistic molecular dynamics (MD) simulations are routinely used to calculate free energies of small molecules binding to proteins, crossing lipid membranes, and solvation but are computationally expensive. Machine learning (ML) and empirical methods are also used throughout drug discovery but rely on experimental data, limiting the domain of applicability. We present atomistic MD simulations calculating 15,000 small molecule free energies of transfer from water to cyclohexane. This large data set is used to train ML models that predict the free energies of transfer. We show that a spatial graph neural network model achieves the highest accuracy, followed closely by a 3D-convolutional neural network, and shallow learning based on the chemical fingerprint is significantly less accurate. A mean absolute error of ∼4 kJ/mol compared to the MD calculations was achieved for our best ML model. We also show that including data from the MD simulation improves the predictions, tests the transferability of each model to a diverse set of molecules, and show multitask learning improves the predictions. This work provides insight into the hydrophobicity of small molecules and ML cheminformatics modeling, and our data set will be useful for designing and testing future ML cheminformatics methods.

Molecular Mechanism for Gramicidin Dimerization and Dissociation in Bilayers of Different Thickness.

Membrane protein functions can be altered by subtle changes in the host lipid bilayer physical properties. Gramicidin channels have emerged as a powerful system for elucidating the underlying mechanisms of membrane protein function regulation through changes in bilayer properties, which are reflected in the thermodynamic equilibrium distribution between nonconducting gramicidin monomers and conducting bilayer-spanning dimers. To improve our understanding of how subtle changes in bilayer thickness alter the gramicidin monomer and dimer distributions, we performed extensive atomistic molecular dynamics simulations and fluorescence-quenching experiments on gramicidin A (gA). The free-energy calculations predicted a nonlinear coupling between the bilayer thickness and channel formation. The energetic barrier inhibiting gA channel formation was sharply increased in the thickest bilayer (1,2-dierucoyl-sn-glycero-3-phosphocholine). This prediction was corroborated by experimental results on gramicidin channel activity in bilayers of different thickness. To further explore the mechanism of channel formation, we performed extensive unbiased molecular dynamics simulations, which allowed us to observe spontaneous gA dimer formation in lipid bilayers. The simulations revealed structural rearrangements in the gA subunits and changes in lipid packing, as well as water reorganization, that occur during the dimerization process. Together, the simulations and experiments provide new, to our knowledge, insights into the process and mechanism of gramicidin channel formation, as a prototypical example of the bilayer regulation of membrane protein function.

Phospholipid Chain Interactions with Cholesterol Drive Domain Formation in Lipid Membranes.

Cholesterol is a key component of eukaryotic membranes, but its role in cellular biology in general and in lipid rafts in particular remains controversial. Model membranes are used extensively to determine the phase behavior of ternary mixtures of cholesterol, a saturated lipid, and an unsaturated lipid with liquid-ordered and liquid-disordered phase coexistence. Despite many different experiments that determine lipid-phase diagrams, we lack an understanding of the molecular-level driving forces for liquid phase coexistence in bilayers with cholesterol. Here, we use atomistic molecular dynamics computer simulations to address the driving forces for phase coexistence in ternary lipid mixtures. Domain formation is directly observed in a long-timescale simulation of a mixture of 1,2-distearoyl-sn-glycero-3-phosphocholine, unsaturated 1,2-dilinoleoyl-sn-glycero-3-phosphocholine, and cholesterol. Free-energy calculations for the exchange of the saturated and unsaturated lipids between the ordered and disordered phases give insight into the mixing behavior. We show that a large energetic contribution to domain formation is favorable enthalpic interactions of the saturated lipid in the ordered phase. This favorable energy for forming an ordered, cholesterol-rich phase is opposed by a large unfavorable entropy. Martini coarse-grained simulations capture the unfavorable free energy of mixing but do not reproduce the entropic contribution because of the reduced representation of the phospholipid tails. Phospholipid tails and their degree of unsaturation are key energetic contributors to lipid phase separation.

The importance of membrane defects-lessons from simulations.

The defects and pores within lipid membranes are scientifically interesting and have a number of biological applications. Although lipid bilayers are extremely thin hydrophobic barriers, just ∼3 nm thick, they include diverse chemistry and have complex structures. Bilayers are soft and dynamic, and as a result, they can bend and deform in response to different stimuli by means of structural changes in their component lipids. Though defects occur within these structures, their transience and small size have made it difficult to characterize them. However, with recent advances in computer power and computational modeling techniques, researchers can now use simulations as a powerful tool to probe the mechanism and energies of defect and pore formation in a number of situations. In this Account, we present results from our detailed molecular dynamics computer simulations of hydrophilic pores and related defects in lipid bilayers at an atomistic level. Electroporation can be used to increase the permeability of cellular membranes, with potential therapeutic applications. Atomistic simulations of electroporation have illustrated the molecular details of this process, including the importance of water dipole interactions at the water-membrane interface. Characterization of the lipid-protein interactions provides an important tool for understanding transmembrane protein structure and thermodynamic stability. Atomistic simulations give a detailed picture of the free energies of model peptides and side chains in lipid membranes; the energetic cost of defect formation strongly influences the energies of interactions between lipids and polar and charged residues. Many antimicrobial peptides form hydrophilic pores in lipid membranes, killing bacteria or cancer cells. On the basis of simulation data, at least some of these peptides form defects and pores near the center of the bilayer, with a common disordered structure where hydrated headgroups form an approximately toroidal shape. The localization and trafficking of lipids supports general membrane structure and a number of important signaling cascades, such as those involving ceramide, diacylglycerol, and cholesterol. Atomistic simulations have determined the rates and free energies of lipid flip-flop. During the flip-flop of most phosphatidylcholine lipids, a hydrophilic pore forms when the headgroup moves near the center of the bilayer. Simulations have provided novel insight into many features of defects and pores in lipid membranes. Simulation data from very different systems and models show how water penetration and defect formation can determine the free energies of many membrane processes. Bilayers can deform and allow transient defects and pores when exposed to a diverse range of stimuli. Future work will explore many aspects of membrane defects with increased resolution and scope, including the study of more complex lipid mixtures, membrane domains, and large-scale membrane remodeling. Such studies will examine processes including vesicle budding and fusion, non-bilayer lipid phases, and interactions between lipid bilayers and other biomolecules. Simulations provide information that complements experimental studies, allowing microscopic insight into experimental observations and suggesting novel hypotheses and experiments. These studies should enable a deeper understanding of the role of lipid bilayers in cellular biology and support the development of future lipid-based biotechnology.

Folding and insertion thermodynamics of the transmembrane WALP peptide.

The anchor of most integral membrane proteins consists of one or several helices spanning the lipid bilayer. The WALP peptide, GWW(LA)n (L)WWA, is a common model helix to study the fundamentals of protein insertion and folding, as well as helix-helix association in the membrane. Its structural properties have been illuminated in a large number of experimental and simulation studies. In this combined coarse-grained and atomistic simulation study, we probe the thermodynamics of a single WALP peptide, focusing on both the insertion across the water-membrane interface, as well as folding in both water and a membrane. The potential of mean force characterizing the peptide's insertion into the membrane shows qualitatively similar behavior across peptides and three force fields. However, the Martini force field exhibits a pronounced secondary minimum for an adsorbed interfacial state, which may even become the global minimum-in contrast to both atomistic simulations and the alternative PLUM force field. Even though the two coarse-grained models reproduce the free energy of insertion of individual amino acids side chains, they both underestimate its corresponding value for the full peptide (as compared with atomistic simulations), hinting at cooperative physics beyond the residue level. Folding of WALP in the two environments indicates the helix as the most stable structure, though with different relative stabilities and chain-length dependence.

Atomistic simulations of pore formation and closure in lipid bilayers.

Cellular membranes separate distinct aqueous compartments, but can be breached by transient hydrophilic pores. A large energetic cost prevents pore formation, which is largely dependent on the composition and structure of the lipid bilayer. The softness of bilayers and the disordered structure of pores make their characterization difficult. We use molecular-dynamics simulations with atomistic detail to study the thermodynamics, kinetics, and mechanism of pore formation and closure in DLPC, DMPC, and DPPC bilayers, with pore formation free energies of 17, 45, and 78 kJ/mol, respectively. By using atomistic computer simulations, we are able to determine not only the free energy for pore formation, but also the enthalpy and entropy, which yields what is believed to be significant new insights in the molecular driving forces behind membrane defects. The free energy cost for pore formation is due to a large unfavorable entropic contribution and a favorable change in enthalpy. Changes in hydrogen bonding patterns occur, with increased lipid-water interactions, and fewer water-water hydrogen bonds, but the total number of overall hydrogen bonds is constant. Equilibrium pore formation is directly observed in the thin DLPC lipid bilayer. Multiple long timescale simulations of pore closure are used to predict pore lifetimes. Our results are important for biological applications, including the activity of antimicrobial peptides and a better understanding of membrane protein folding, and improve our understanding of the fundamental physicochemical nature of membranes.

Computer simulations of lipid membrane domains.

There is great diversity in the composition and structure of biological lipid membranes. We are beginning to appreciate the crucial role of lipids in many cellular processes, and characterize some of the lateral structures within membranes that could play a role in the activity of lipids. Simulations probe molecular level interactions between single molecules, which provide complementary information to experiments. Lipid membrane simulations have reached an exciting point, where the time and length scales of our simulations are approaching experimental resolutions and can be used to interpret experiments on increasingly complex model membranes. The focus of this review is on recent molecular simulations of domain formation in lipid bilayers. We highlight a number of recent examples where simulations are used in collaboration with experiments. We review recent simulation studies on lipid-lipid interactions related to domain formation and on lipid-protein interactions relevant for lipid raft function.

Oleic acid phase behavior from molecular dynamics simulations.

Fatty acid aggregation is important for a number of diverse applications: from origins of life research to industrial applications to health and disease. Experiments have characterized the phase behavior of oleic acid mixtures, but the molecular details are complex and hard to probe with many experiments. Coarse-grained molecular dynamics computer simulations and free energy calculations are used to model oleic acid aggregation. From random dispersions, we observe the aggregation of oleic acid monomers into micelles, vesicles, and oil phases, depending on the protonation state of the oleic acid head groups. Worm-like micelles are observed when all the oleic acid molecules are deprotonated and negatively charged. Vesicles form spontaneously if significant amounts of both neutral and negative oleic acid are present. Oil phases form when all the fatty acids are protonated and neutral. This behavior qualitatively matches experimental observations of oleic acid aggregation. To explain the observed phase behavior, we use umbrella sampling free energy calculations to determine the stability of individual monomers in aggregates compared to water. We find that both neutral and negative oleic acid molecules prefer larger aggregates, but neutral monomers prefer negatively charged aggregates and negative monomers prefer neutral aggregates. Both neutral and negative monomers are most stable in a DOPC bilayer, with implications on fatty acid adsorption and cellular membrane evolution. Although the CG model qualitatively reproduces oleic acid phase behavior, we show that an updated polarizable water model is needed to more accurately predict the shift in pKa for oleic acid in model bilayers.

Conical lipids in flat bilayers induce packing defects similar to that induced by positive curvature.

In biological membranes, changes in lipid composition or mechanical deformations produce defects in the geometrical arrangement of lipids, thus allowing the adsorption of certain peripheral proteins. Here, we perform molecular dynamics simulations on bilayers containing a cylindrical lipid (PC) and a conical lipid (DOG). Profiles of atomic density and lateral pressure across the bilayer show differences in the acyl chain region due to deeper partitioning of DOG compared to PC. However, such analyses are less informative for the interfacial region where peripheral proteins adsorb. To circumvent this limitation, we develop, to our knowledge, a new method of membrane surface analysis. This method allows the identification of chemical defects, where hydrocarbon chains are accessible to the solvent, and geometrical defects, i.e., voids deeper than the glycerol backbone. The size and number of both types of defects increase with the number of monounsaturated acyl chains in PC and with the introduction of DOG, although the defects do not colocalize with the conical lipid. Interestingly, the size and probability of the defects promoted by DOG resemble those induced by positive curvature, thus explaining why conical lipids and positive curvature can both drive the adsorption of peripheral proteins that use hydrophobic residues as membrane anchors.

Advances in Computational Approaches for Estimating Passive Permeability in Drug Discovery.

Passive permeation of cellular membranes is a key feature of many therapeutics. The relevance of passive permeability spans all biological systems as they all employ biomembranes for compartmentalization. A variety of computational techniques are currently utilized and under active development to facilitate the characterization of passive permeability. These methods include lipophilicity relations, molecular dynamics simulations, and machine learning, which vary in accuracy, complexity, and computational cost. This review briefly introduces the underlying theories, such as the prominent inhomogeneous solubility diffusion model, and covers a number of recent applications. Various machine-learning applications, which have demonstrated good potential for high-volume, data-driven permeability predictions, are also discussed. Due to the confluence of novel computational methods and next-generation exascale computers, we anticipate an exciting future for computationally driven permeability predictions.

Molecular simulation of rapid translocation of cholesterol, diacylglycerol, and ceramide in model raft and nonraft membranes.

The translocation of lipids across membranes (flip-flop) is an important biological process. Slow exchange on a physiological timescale allows the creation of asymmetric distributions of lipids across cellular membranes. The location of lipids and their rate of exchange have important biological consequences, especially for lipids involved in cellular signaling. We investigated the translocation of cholesterol, ceramide, and diacylglycerol in two model bilayers using molecular dynamics simulations. We estimate half times for flip-flop for cholesterol, diacylglycerol, and ceramide of 20 µs, 30 µs, and 10 ms in a POPC bilayer, compared with approximately 30 min, 30 ms, and 30 s in a model raft bilayer (1:1:1 PSM, POPC, and cholesterol). Cholesterol has a large (54 kJ/mol) free energy of exchange between the POPC and raft bilayer, and therefore, it strongly prefers the more ordered and rigid raft bilayer over the more liquid POPC bilayer. Ceramide and diacylglycerol have relatively small free energies of exchange, suggesting nearly equal preference for both bilayers. This unexpected result may have implications for ceramide and diacylglycerol signaling and membrane localization.

Structural arrangement of the transmission interface in the antigen ABC transport complex TAP.

The transporter associated with antigen processing (TAP) represents a focal point in the immune recognition of virally or malignantly transformed cells by translocating proteasomal degradation products into the endoplasmic reticulum-lumen for loading of MHC class I molecules. Based on a number of experimental data and the homology to the bacterial ABC exporter Sav1866, we constructed a 3D structural model of the core TAP complex and used it to examine the interface between the transmembrane and nucleotide-binding domains (NBD) by cysteine-scanning and cross-linking approaches. Herein, we demonstrate the functional importance of the newly identified X-loop in the NBD in coupling substrate binding to downstream events in the transport cycle. We further verified domain swapping in a heterodimeric ABC half-transporter complex by cysteine cross-linking. Strikingly, either substrate binding or translocation can be blocked by cross-linking the X-loop to coupling helix 2 or 1, respectively. These results resolve the structural arrangement of the transmission interface and point to different functions of the cytosolic loops and coupling helices in substrate binding, signaling, and transport.

Transfer of arginine into lipid bilayers is nonadditive.

Computer simulations suggest that the translocation of arginine through the hydrocarbon core of a lipid membrane proceeds by the formation of a water-filled defect that keeps the arginine molecule hydrated even at the center of the bilayer. We show here that adding additional arginine molecules into one of these water defects causes only a small change in free energy. The barrier for transferring multiple arginines through the membrane is approximately the same as for a single arginine and may even be lower depending on the exact geometry of the system. We discuss these results in the context of arginine-rich peptides such as antimicrobial and cell-penetrating peptides.

Curvature Energetics Determined by Alchemical Simulation on Four Topologically Distinct Lipid Phases.

The relative curvature energetics of two lipids are tested using thermodynamic integration (TI) on four topologically distinct lipid phases. Simulations use TI to switch between choline headgroup lipids (POPC; that prefers to be flat) and ethanolamine headgroup lipids (POPE; that prefer, for example, the inner monolayer of vesicles). The thermodynamical moving of the lipids between planar, inverse hexagonal (HII), cubic (QII; Pn3m space group), and vesicle topologies reveals differences in material parameters that were previously challenging to access. The methodology allows for predictions of two important lipid material properties: the difference in POPC/POPE monolayer intrinsic curvature (ΔJ0) and the difference in POPC/POPE monolayer Gaussian curvature modulus (Δκ̅m), both of which are connected to the energetics of topological variation. Analysis of the TI data indicates that, consistent with previous experiment and simulation, the J0 of POPE is more negative than POPC (ΔJ0 = -0.018 ± 0.001 Å-1). The theoretical framework extracts significant differences in κ̅m of which POPE is less negative than POPC by 2.0 to 4.0 kcal/mol. The range of these values is determined by considering subsets of the simulations, and disagreement between these subsets suggests separate mechanical parameters at very high curvature. Finally, the fit of the TI data to the model indicates that the position of the pivotal plane of curvature is not constant across topologies at high curvature. Overall, the results offer insights into lipid material properties, the limits of a single HC model, and how to test them using simulation.

Atomistic Characterization of Gramicidin Channel Formation.

We investigated gramicidin A (gA) subunit dimerization in lipid bilayers using microsecond-long replica-exchange umbrella sampling simulations, millisecond-long unbiased molecular dynamics simulations, and machine learning. Our simulations led to a dimer structure that is indistinguishable from the experimentally determined gA channel structures, with the two gA subunits joined by six hydrogen bonds (6HB). The simulations also uncovered two additional dimer structures, with different gA-gA stacking orientations that were stabilized by four or two hydrogen bonds (4HB or 2HB). When examining the temporal evolution of the dimerization, we found that two bilayer-inserted gA subunits can form the 6HB dimer directly, with no discernible intermediate states, as well as through paths that involve the 2HB and 4HB dimers.

Discovery of Small-Molecule Inhibitors of SARS-CoV-2 Proteins Using a Computational and Experimental Pipeline.

A rapid response is necessary to contain emergent biological outbreaks before they can become pandemics. The novel coronavirus (SARS-CoV-2) that causes COVID-19 was first reported in December of 2019 in Wuhan, China and reached most corners of the globe in less than two months. In just over a year since the initial infections, COVID-19 infected almost 100 million people worldwide. Although similar to SARS-CoV and MERS-CoV, SARS-CoV-2 has resisted treatments that are effective against other coronaviruses. Crystal structures of two SARS-CoV-2 proteins, spike protein and main protease, have been reported and can serve as targets for studies in neutralizing this threat. We have employed molecular docking, molecular dynamics simulations, and machine learning to identify from a library of 26 million molecules possible candidate compounds that may attenuate or neutralize the effects of this virus. The viability of selected candidate compounds against SARS-CoV-2 was determined experimentally by biolayer interferometry and FRET-based activity protein assays along with virus-based assays. In the pseudovirus assay, imatinib and lapatinib had IC50 values below 10 µM, while candesartan cilexetil had an IC50 value of approximately 67 µM against Mpro in a FRET-based activity assay. Comparatively, candesartan cilexetil had the highest selectivity index of all compounds tested as its half-maximal cytotoxicity concentration 50 (CC50) value was the only one greater than the limit of the assay (>100 µM).

Molecular simulations of lipid flip-flop in the presence of model transmembrane helices.

The transport of lipids between membrane leaflets, also known as flip-flop, is a key process in regulating the lipid composition of biological membranes. It is also important for the growth of biogenic membranes that are the site for lipid synthesis. It has been shown that the mere presence of transmembrane alpha-helical peptides or proteins enhances the rate lipid flip-flop [Kol et al. (2001) Biochemistry 40, 10500-10506]. Using computational models of natural phospholipids with different headgroups, we calculated the free energy profiles for transferring single phospholipids from bulk water to the center of a dioleylphosphatidylcholine (DOPC) bilayer in the presence of transmembrane helices. The free energy barrier for phosphatidylethanolamine (PE) and phosphatidylglycerol (PG) flip-flop decreased by a few kilojoules per mole when a WALP23 or KALP23 peptide was present in the membrane, while the barrier for PC was not affected. We observed large bilayer deformations during lipid flip-flop when the hydrophilic headgroup is in the hydrophobic interior of the bilayer. The presence of KALP23 or WALP23 decreased the size and stability of these defects, suggesting integral membrane proteins affect the mechanism of flip-flop. There was a large decrease in the free energy of desorption for PE and PG when transmembrane peptides were present. This suggests specific PE and PG interactions with the peptide have a large affect on their stability in the membrane, with implications on cellular lipid and protein trafficking.