GPU-specific algorithms for improved solute sampling in grand canonical Monte Carlo simulations.

The Grand Canonical Monte Carlo (GCMC) ensemble defined by the excess chemical potential, µex , volume, and temperature, in the context of molecular simulations allows for variations in the number of particles in the system. In practice, GCMC simulations have been widely applied for the sampling of rare gasses and water, but limited in the context of larger molecules. To overcome this limitation, the oscillating µex GCMC method was introduced and shown to be of utility for sampling small solutes, such as formamide, propane, and benzene, as well as for ionic species such as monocations, acetate, and methylammonium. However, the acceptance of GCMC insertions is low, and the method is computationally demanding. In the present study, we improved the sampling efficiency of the GCMC method using known cavity-bias and configurational-bias algorithms in the context of GPU architecture. Specifically, for GCMC simulations of aqueous solution systems, the configurational-bias algorithm was extended by applying system partitioning in conjunction with a random interval extraction algorithm, thereby improving the efficiency in a highly parallel computing environment. The method is parallelized on the GPU using CUDA and OpenCL, allowing for the code to run on both Nvidia and AMD GPUs, respectively. Notably, the method is particularly well suited for GPU computing as the large number of threads allows for simultaneous sampling of a large number of configurations during insertion attempts without additional computational overhead. In addition, the partitioning scheme allows for simultaneous insertion attempts for large systems, offering considerable efficiency. Calculations on the BK Channel, a transporter, including a lipid bilayer with over 760,000 atoms, show a speed up of ~53-fold through the use of system partitioning. The improved algorithm is then combined with an enhanced µex oscillation protocol and shown to be of utility in the context of the site-identification by ligand competitive saturation (SILCS) co-solvent sampling approach as illustrated through application to the protein CDK2.

Uncovering the folding mechanism of pertactin: A comparative study of isolated and vectorial folding.

Autotransporters are a large family of virulence factors found in Gram-negative bacteria that play important roles in their pathogenesis. The passenger domain of autotransporters is almost always composed of a large ß-helix, with only a small portion of it being relevant to its virulence function. This has led to the hypothesis that the folding of the ß-helical structure aids the secretion of the passenger domain across the Gram-negative outer membrane. In this study, we used molecular dynamics simulations and enhanced sampling methods to investigate the stability and folding of the passenger domain of pertactin, an autotransporter from Bordetella pertussis. Specifically, we employed steered molecular dynamics to simulate the unfolding of the entire passenger domain as well as self-learning adaptive umbrella sampling to compare the energetics of folding rungs of the ß-helix independently ("isolated folding") versus folding rungs on top of a previously folded rung ("vectorial folding"). Our results showed that vectorial folding is highly favorable compared with isolated folding; moreover, our simulations showed that the C-terminal rung of the ß-helix is the most resistant to unfolding, in agreement with previous studies that found the C-terminal half of the passenger domain to be more stable than the N-terminal one. Overall, this study provides new insights into the folding process of an autotransporter passenger domain and its potential role in secretion across the outer membrane.

Considerations for Use of Pharmacodynamic Biomarkers to Support Biosimilar Development - (III) A Randomized Trial with Interferon Beta-1a Products.

The US Food and Drug Administration (FDA) has taken steps to bring efficiency to the development of biosimilars, including establishing guidance for the use of pharmacokinetic and pharmacodynamic (PD) similarity study data without a comparative clinical study with efficacy end point(s). To better understand the potential role for PD biomarkers in biosimilar development and inform best practices for biomarker selection and analysis, we conducted a randomized, double-blinded, placebo-controlled, single-dose, parallel-arm clinical study in healthy participants. Eighty-four healthy participants (n = 12 per dose arm) received either placebo or one of three doses of either interferon ß-1a (7.5-30 µg) or pegylated interferon ß-1a (31.25-125 µg) to evaluate the maximum change from baseline and the baseline-adjusted area under the effect curve for the biomarkers neopterin in serum and myxovirus resistance protein 1 in blood. Both PD biomarkers increased following product administration with clear separation from baseline (neopterin: 3.4-fold and 3.9-fold increase for interferon ß-1a and pegylated interferon ß-1a, respectively; myxovirus resistance protein 1: 19.0-fold and 47.2-fold increase for interferon ß-1a and pegylated interferon ß-1a, respectively). The dose-response curves support that therapeutic doses were adequately sensitive to detect differences in both PD biomarkers for consideration in a PD similarity study design. Because baseline levels of both biomarkers are low compared with on-treatment values, there was little difference in using PD measures adjusted to baseline compared with the results without baseline adjustment. This study illustrates potential methodologies for evaluating PD biomarkers and an approach to address information gaps when limited information is publicly available for one or more PD biomarkers.

SILCS-RNA: Toward a Structure-Based Drug Design Approach for Targeting RNAs with Small Molecules.

RNA molecules can act as potential drug targets in different diseases, as their dysregulated expression or misfolding can alter various cellular processes. Noncoding RNAs account for ∼70% of the human genome, and these molecules can have complex tertiary structures that present a great opportunity for targeting by small molecules. In the present study, the site identification by ligand competitive saturation (SILCS) computational approach is extended to target RNA, termed SILCS-RNA. Extensions to the method include an enhanced oscillating excess chemical potential protocol for the grand canonical Monte Carlo calculations and individual simulations of the neutral and charged solutes from which the SILCS functional group affinity maps (FragMaps) are calculated for subsequent binding site identification and docking calculations. The method is developed and evaluated against seven RNA targets and their reported small molecule ligands. SILCS-RNA provides a detailed characterization of the functional group affinity pattern in the small molecule binding sites, recapitulating the types of functional groups present in the ligands. The developed method is also shown to be useful for identification of new potential binding sites and identifying ligand moieties that contribute to binding, granular information that can facilitate ligand design. However, limitations in the method are evident including the ability to map the regions of binding sites occupied by ligand phosphate moieties and to fully account for the wide range of conformational heterogeneity in RNA associated with binding of different small molecules, emphasizing inherent challenges associated with applying computer-aided drug design methods to RNA. While limitations are present, the current study indicates how the SILCS-RNA approach may enhance drug discovery efforts targeting RNAs with small molecules.

Application of Site-Identification by Ligand Competitive Saturation in Computer-Aided Drug Design.

Site Identification by Ligand Competitive Saturation (SILCS) is a molecular simulation approach that uses diverse small solutes in aqueous solution to obtain functional group affinity patterns of a protein or other macromolecule. This involves employing a combined Grand Canonical Monte Carlo (GCMC)-molecular dynamics (MD) method to sample the full 3D space of the protein, including deep binding pockets and interior cavities from which functional group free energy maps (FragMaps) are obtained. The information content in the maps, which include contributions from protein flexibilty and both protein and functional group desolvation contributions, can be used in many aspects of the drug discovery process. These include identification of novel ligand binding pockets, including allosteric sites, pharmacophore modeling, prediction of relative protein-ligand binding affinities for database screening and lead optimization efforts, evaluation of protein-protein interactions as well as in the formulation of biologics-based drugs including monoclonal antibodies. The present article summarizes the various tools developed in the context of the SILCS methodology and their utility in computer-aided drug design (CADD) applications, showing how the SILCS toolset can improve the drug-development process on a number of fronts with respect to both accuracy and throughput representing a new avenue of CADD applications.

Lpp positions peptidoglycan at the AcrA-TolC interface in the AcrAB-TolC multidrug efflux pump.

The multidrug efflux pumps of Gram-negative bacteria are a class of complexes that span the periplasm, coupling both the inner and outer membranes to expel toxic molecules. The best-characterized example of these tripartite pumps is the AcrAB-TolC complex of Escherichia coli. However, how the complex interacts with the peptidoglycan (PG) cell wall, which is anchored to the outer membrane (OM) by Braun's lipoprotein (Lpp), is still largely unknown. In this work, we present molecular dynamics simulations of a complete, atomistic model of the AcrAB-TolC complex with the inner membrane, OM, and PG layers all present. We find that the PG localizes to the junction of AcrA and TolC, in agreement with recent cryo-tomography data. Free-energy calculations reveal that the positioning of PG is determined by the length and conformation of multiple Lpp copies anchoring it to the OM. The distance between the PG and OM measured in cryo-electron microscopy images of wild-type E. coli also agrees with the simulation-derived spacing. Sequence analysis of AcrA suggests a conserved role for interactions with PG in the assembly and stabilization of efflux pumps, one that may extend to other trans-envelope complexes as well.

Rapid and accurate estimation of protein-ligand relative binding affinities using site-identification by ligand competitive saturation.

Predicting relative protein-ligand binding affinities is a central pillar of lead optimization efforts in structure-based drug design. The site identification by ligand competitive saturation (SILCS) methodology is based on functional group affinity patterns in the form of free energy maps that may be used to compute protein-ligand binding poses and affinities. Presented are results obtained from the SILCS methodology for a set of eight target proteins as reported originally in Wang et al. (J. Am. Chem. Soc., 2015, 137, 2695-2703) using free energy perturbation (FEP) methods in conjunction with enhanced sampling and cycle closure corrections. These eight targets have been subsequently studied by many other authors to compare the efficacy of their method while comparing with the outcomes of Wang et al. In this work, we present results for a total of 407 ligands on the eight targets and include specific analysis on the subset of 199 ligands considered previously. Using the SILCS methodology we can achieve an average accuracy of up to 77% and 74% when considering the eight targets with their 199 and 407 ligands, respectively, for rank-ordering ligand affinities as calculated by the percent correct metric. This accuracy increases to 82% and 80%, respectively, when the SILCS atomic free energy contributions are optimized using a Bayesian Markov-chain Monte Carlo approach. We also report other metrics including Pearson's correlation coefficient, Pearlman's predictive index, mean unsigned error, and root mean square error for both sets of ligands. The results obtained for the 199 ligands are compared with the outcomes of Wang et al. and other published works. Overall, the SILCS methodology yields similar or better-quality predictions without a priori need for known ligand orientations in terms of the different metrics when compared to current FEP approaches with significant computational savings while additionally offering quantitative estimates of individual atomic contributions to binding free energies. These results further validate the SILCS methodology as an accurate, computationally efficient tool to support lead optimization and drug discovery.

A Minimal Membrane Metal Transport System: Dynamics and Energetics of mer Proteins.

The mer operon in bacteria encodes a set of proteins and enzymes that impart resistance to environmental mercury toxicity by importing Hg2+ and reducing it to volatile Hg(0). Because the reduction occurs in the cytoplasm, mercuric ions must first be transported across the cytoplasmic membrane by one of a few known transporters. MerF is the smallest of these, containing only two transmembrane helices and two pairs of vicinal cysteines that coordinate mercuric ions. In this work, we use molecular dynamics simulations to characterize the dynamics of MerF in its apo and Hg2+ -bound states. We find that the apo state positions one of the cysteine pairs closer to the periplasmic side of the membrane, while in the bound state the same pair approaches the cytoplasmic side. This finding is consistent with the functional requirement of accepting Hg2+ from the periplasmic space, sequestering it on acceptance, and transferring it to the cytoplasm. Conformational changes in the TM helices facilitate the functional interaction of the two cysteine pairs. Free-energy calculations provide a barrier of 16 kcal/mol for the association of the periplasmic Hg2+ -bound protein MerP with MerF and 7 kcal/mol for the subsequent association of MerF's two cysteine pairs. Despite the significant conformational changes required to move the binding site across the membrane, coarse-grained simulations of multiple copies of MerF support the expectation that it functions as a monomer. Our results demonstrate how conformational changes and binding thermodynamics could lead to such a small membrane protein acting as an ion transporter. Published 2019. This article is a U.S. Government work and is in the public domain in the USA.

Conformational Dynamics of AcrA Govern Multidrug Efflux Pump Assembly.

Multidrug efflux pumps of pathogenic, Gram-negative bacteria comprise an innate resistance mechanism and are key contributors to the emerging global pandemic of antibiotic resistance. Several increasingly detailed cryo-electron microscopy maps have been resolved of an entire efflux pump complex, AcrAB-TolC, resulting in atomistic structural models. Using a recent model, we have carried out nearly 40 µs of molecular dynamics simulations to study one of the key components of the protein complex AcrA, the membrane fusion protein that connects the inner-membrane-bound AcrB to the outer-membrane-bound TolC. We determined a three-dimensional potential of mean force (PMF) for AcrA, which displays two main conformational basins representing assembly competent and incompetent states. Corresponding experiments show that stabilizing mutations at an interdomain interface shift the dynamic equilibrium between these states to the incompetent one, disrupting pump assembly and function and resensitizing bacteria to existing antibiotics. The modulation of AcrA dynamics through pharmacological intervention therefore presents a promising route for the development of new antibiotics.

Stable calcium-free myocilin olfactomedin domain variants reveal challenges in differentiating between benign and glaucoma-causing mutations.

Nonsynonymous gene mutations can be beneficial, neutral, or detrimental to the stability, structure, and biological function of the encoded protein, but the effects of these mutations are often not readily predictable. For example, the ß-propeller olfactomedin domain of myocilin (mOLF) exhibits a complex interrelationship among structure(s), stability, and aggregation. Numerous mutations within mOLF are linked to glaucoma; the resulting variants are less stable, aggregation-prone, and sequestered intracellularly, causing cytotoxicity. Here, we report the first stable mOLF variants carrying substitutions in the calcium-binding site that exhibit solution characteristics indistinguishable from those of glaucoma variants. Crystal structures of these stable variants at 1.8-2.0-Å resolution revealed features that we could not predict by molecular dynamics simulations, including loss of loop structure, helix unwinding, and a blade shift. Double mutants that combined a stabilizing substitution and a selected glaucoma-causing single-point mutant rescued in vitro folding and stability defects. In the context of full-length myocilin, secretion of stable single variants was indistinguishable from that of the WT protein, and the double mutants were secreted to varying extents. In summary, our finding that mOLF can tolerate particular substitutions that render the protein stable despite a conformational switch emphasizes the complexities in differentiating between benign and glaucoma-causing variants and provides new insight into the possible biological function of myocilin.

Accelerating Membrane Simulations with Hydrogen Mass Repartitioning.

The time step of atomistic molecular dynamics (MD) simulations is determined by the fastest motions in the system and is typically limited to 2 fs. An increasingly popular approach is to increase the mass of the hydrogen atoms to ∼3 amu and decrease the mass of the parent atom by an equivalent amount. This approach, known as hydrogen-mass repartitioning (HMR), permits time steps up to 4 fs with reasonable simulation stability. While HMR has been applied in many published studies to date, it has not been extensively tested for membrane-containing systems. Here, we compare the results of simulations of a variety of membranes and membrane-protein systems run using a 2 fs time step and a 4 fs time step with HMR. For pure membrane systems, we find almost no difference in structural properties, such as area-per-lipid, electron density profiles, and order parameters, although there are differences in kinetic properties such as the diffusion constant. Conductance through a porin in an applied field, partitioning of a small peptide, hydrogen-bond dynamics, and membrane mixing show very little dependence on HMR and the time step. We also tested a 9 Å cutoff as compared to the standard CHARMM cutoff of 12 Å, finding significant deviations in many properties tested. We conclude that HMR is a valid approach for membrane systems, but a 9 Å cutoff is not.

Identification of Binding Sites for Efflux Pump Inhibitors of the AcrAB-TolC Component AcrA.

The overexpression of multidrug efflux pumps is an important mechanism of clinical resistance in Gram-negative bacteria. Recently, four small molecules were discovered that inhibit efflux in Escherichia coli and interact with the AcrAB-TolC efflux pump component AcrA. However, the binding site(s) for these molecules was not determined. Here, we combine ensemble docking and molecular dynamics simulations with tryptophan fluorescence spectroscopy, site-directed mutagenesis, and antibiotic susceptibility assays to probe binding sites and effects of binding of these molecules. We conclude that clorobiocin and SLU-258 likely bind at a site located between the lipoyl and ß-barrel domains of AcrA.

Transmembrane but not soluble helices fold inside the ribosome tunnel.

Integral membrane proteins are assembled into the ER membrane via a continuous ribosome-translocon channel. The hydrophobicity and thickness of the core of the membrane bilayer leads to the expectation that transmembrane (TM) segments minimize the cost of harbouring polar polypeptide backbones by adopting a regular pattern of hydrogen bonds to form α-helices before integration. Co-translational folding of nascent chains into an α-helical conformation in the ribosomal tunnel has been demonstrated previously, but the features governing this folding are not well understood. In particular, little is known about what features influence the propensity to acquire α-helical structure in the ribosome. Using in vitro translation of truncated nascent chains trapped within the ribosome tunnel and molecular dynamics simulations, we show that folding in the ribosome is attained for TM helices but not for soluble helices, presumably facilitating SRP (signal recognition particle) recognition and/or a favourable conformation for membrane integration upon translocon entry.

Folding free energy landscapes of ß-sheets with non-polarizable and polarizable CHARMM force fields.

Molecular dynamics (MD) simulations of peptides and proteins offer atomic-level detail into many biological processes, although the degree of insight depends on the accuracy of the force fields used to represent them. Protein folding is a key example in which the accurate reproduction of folded-state conformations of proteins and kinetics of the folding processes in simulation is a longstanding goal. Although there have been a number of recent successes, challenges remain in capturing the full complexity of folding for even secondary-structure elements. In the present work, we have used all-atom MD simulations to study the folding properties of one such element, the C-terminal ß-hairpin of the B1 domain of streptococcal protein G (GB1). Using replica-exchange umbrella sampling simulations, we examined the folding free energy of two fixed-charge CHARMM force fields, CHARMM36 and CHARMM22*, as well as a polarizable force field, the CHARMM Drude-2013 model, which has previously been shown to improve the folding properties of α-helical peptides. The CHARMM22* and Drude-2013 models are in rough agreement with experimental studies of GB1 folding, while CHARMM36 overstabilizes the ß-hairpin. Additional free-energy calculations show that small adjustments to the atomic polarizabilities in the Drude-2013 model can improve both the backbone solubility and folding properties of GB1 without significantly affecting the model's ability to properly fold α-helices. We also identify a non-native salt bridge in the ß-turn region that overstabilizes the ß-hairpin in the C36 model. Finally, we demonstrate that tryptophan fluorescence is insufficient for capturing the full ß-hairpin folding pathway.

Correction to: Computed Free Energies of Peptide Insertion into Bilayers are Independent of Computational Method.

The original version of the article unfortunately contained an error in NIH support grant number RO1-GM74639 in the Acknowledgements section. The correct grant number is RO1-GM74637. This has been corrected with this erratum.

Computed Free Energies of Peptide Insertion into Bilayers are Independent of Computational Method.

We show that the free energy of inserting hydrophobic peptides into lipid bilayer membranes from surface-aligned to transmembrane inserted states can be reliably calculated using atomistic models. We use two entirely different computational methods: high temperature spontaneous peptide insertion calculations as well as umbrella sampling potential-of-mean-force (PMF) calculations, both yielding the same energetic profiles. The insertion free energies were calculated using two different protein and lipid force fields (OPLS protein/united-atom lipids and CHARMM36 protein/all-atom lipids) and found to be independent of the simulation parameters. In addition, the free energy of insertion is found to be independent of temperature for both force fields. However, we find major difference in the partitioning kinetics between OPLS and CHARMM36, likely due to the difference in roughness of the underlying free energy surfaces. Our results demonstrate not only a reliable method to calculate insertion free energies for peptides, but also represent a rare case where equilibrium simulations and PMF calculations can be directly compared.

Structure and Misfolding of the Flexible Tripartite Coiled-Coil Domain of Glaucoma-Associated Myocilin.

Glaucoma-associated myocilin is a member of the olfactomedins, a protein family involved in neuronal development and human diseases. Molecular studies of the myocilin N-terminal coiled coil demonstrate a unique tripartite architecture: a Y-shaped parallel dimer-of-dimers with distinct tetramer and dimer regions. The structure of the dimeric C-terminal 7-heptad repeats elucidates an unexpected repeat pattern involving inter-strand stabilization by oppositely charged residues. Molecular dynamics simulations reveal an alternate accessible conformation in which the terminal inter-strand disulfide limits the extent of unfolding and results in a kinked configuration. By inference, full-length myocilin is also branched, with two pairs of C-terminal olfactomedin domains. Selected variants within the N-terminal region alter the apparent quaternary structure of myocilin but do so without compromising stability or causing aggregation. In addition to increasing our structural knowledge of naturally occurring extracellular coiled coils and biomedically important olfactomedins, this work broadens the scope of protein misfolding in the pathogenesis of myocilin-associated glaucoma.

Redox-Driven Conformational Dynamics in a Photosystem-II-Inspired ß-Hairpin Maquette Determined through Spectroscopy and Simulation.

Tyrosine-based radical transfer plays an important role in photosynthesis, respiration, and DNA synthesis. Radical transfer can occur either by electron transfer (ET) or proton coupled electron transfer (PCET), depending on the pH. Reversible conformational changes in the surrounding protein matrix may control reactivity of radical intermediates. De novo designed Peptide A is a synthetic 18 amino-acid ß-hairpin, which contains a single tyrosine (Y5) and carries out a kinetically significant PCET reaction between Y5 and a cross-strand histidine (H14). In Peptide A, amide II' (CN) changes are observed in the UV resonance Raman (UVRR) spectrum, associated with tyrosine ET and PCET; these bands were attributed previously to a reversible change in secondary structure. Here, we use molecular dynamics simulations to define this conformational change in Peptide A and its H14-to-cyclohexylalanine variant, Peptide C. Three different Y5 charge states, tyrosine (YH), tyrosinate (Y-), and neutral tyrosyl radical (Y·), are considered. The simulations show that Peptide A-YH and A-Y- retain secondary structure and noncovalent interactions, whereas A-Y· is unstable. In contrast, both Peptide C-Y- and Peptide C-Y· are unstable, due to the loss of the Y5-H14 π-π interaction. These simulations are consistent with previous UVRR experimental results on the two ß-hairpins. Furthermore, they demonstrate the ability of simulations using fixed-charge force fields to accurately capture redox-linked conformational dynamics in a ß-strand peptide.

Thermodynamics of Deca-alanine Folding in Water.

The determination of the folding dynamics of polypeptides and proteins is critical in characterizing their functions in biological systems. Numerous computational models and methods have been developed for studying structure formation at the atomic level. Due to its small size and simple structure, deca-alanine is used as a model system in molecular dynamics (MD) simulations. The free energy of unfolding in vacuum has been studied extensively using the end-to-end distance of the peptide as the reaction coordinate. However, few studies have been conducted in the presence of explicit solvent. Previous results show a significant decrease in the free energy of extended conformations in water, but the α-helical state is still notably favored over the extended state. Although sufficient in vacuum, we show that end-to-end distance is incapable of capturing the full complexity of deca-alanine folding in water. Using α-helical content as a second reaction coordinate, we deduce a more descriptive free-energy landscape, one which reveals a second energy minimum in the extended conformations that is of comparable free energy to the α-helical state. Equilibrium simulations demonstrate the relative stability of the extended and α-helical states in water as well as the transition between the two states. This work reveals both the necessity and challenge of determining a proper reaction coordinate to fully characterize a given process.