A Na pump with reduced stoichiometry is up-regulated by brine shrimp in extreme salinities.

Brine shrimp (Artemia) are the only animals to thrive at sodium concentrations above 4 M. Salt excretion is powered by the Na+,K+-ATPase (NKA), a heterodimeric (αß) pump that usually exports 3Na+ in exchange for 2 K+ per hydrolyzed ATP. Artemia express several NKA catalytic α-subunit subtypes. High-salinity adaptation increases abundance of α2KK, an isoform that contains two lysines (Lys308 and Lys758 in transmembrane segments TM4 and TM5, respectively) at positions where canonical NKAs have asparagines (Xenopus α1's Asn333 and Asn785). Using de novo transcriptome assembly and qPCR, we found that Artemia express two salinity-independent canonical α subunits (α1NN and α3NN), as well as two ß variants, in addition to the salinity-controlled α2KK. These ß subunits permitted heterologous expression of the α2KK pump and determination of its CryoEM structure in a closed, ion-free conformation, showing Lys758 residing within the ion-binding cavity. We used electrophysiology to characterize the function of α2KK pumps and compared it to that of Xenopus α1 (and its α2KK-mimicking single- and double-lysine substitutions). The double substitution N333K/N785K confers α2KK-like characteristics to Xenopus α1, and mutant cycle analysis reveals energetic coupling between these two residues, illustrating how α2KK's Lys308 helps to maintain high affinity for external K+ when Lys758 occupies an ion-binding site. By measuring uptake under voltage clamp of the K+-congener 86Rb+, we prove that double-lysine-substituted pumps transport 2Na+ and 1 K+ per catalytic cycle. Our results show how the two lysines contribute to generate a pump with reduced stoichiometry allowing Artemia to maintain steeper Na+ gradients in hypersaline environments.

Mechanism of voltage gating in the voltage-sensing phosphatase Ci-VSP.

Conformational changes in voltage-sensing domains (VSDs) are driven by the transmembrane electric field acting on the protein charges. Yet, the overall energetics and detailed mechanism of this process are not fully understood. Here, we determined free energy and displacement charge landscapes as well as the major conformations visited during a complete functional gating cycle in the isolated VSD of the phosphatase Ci-VSP (Ci-VSD) comprising four transmembrane helices (segments S1 to S4). Molecular dynamics simulations highlight the extent of S4 movements. In addition to the crystallographically determined activated "Up" and resting "Down" states, the simulations predict two Ci-VSD conformations: a deeper resting state ("down-minus") and an extended activated ("up-plus") state. These additional conformations were experimentally probed via systematic cysteine mutagenesis with metal-ion bridges and the engineering of proton conducting mutants at hyperpolarizing voltages. The present results show that these four states are visited sequentially in a stepwise manner during voltage activation, each step translocating one arginine or the equivalent of ∼1 e0 across the membrane electric field, yielding a transfer of ∼3 e0 charges in total for the complete process.

Classifying Protein-Protein Binding Affinity with Free-Energy Calculations and Machine Learning Approaches.

Understanding the intricate phenomenon of neuronal wiring in the brain is of great interest in neuroscience. In the fruit fly, Drosophila melanogaster, the Dpr-DIP interactome has been identified to play an important role in this process. However, experimental data suggest that a merely limited subset of complexes, essentially 57 out of a total of 231, exhibit strong binding affinity. In this work, we sought to identify the residue-level molecular basis underlying the difference in binding affinity using a state-of-the-art methodology consisting of standard binding free-energy calculations with a geometrical route and machine learning (ML) techniques. We determined the binding affinity for two complexes using statistical mechanics simulations, achieving an excellent reproduction of the experimental data. Moreover, we predicted the binding free energy for two additional low-affinity complexes, devoid of experimental estimation, while simultaneously identifying key residues for the binding. Furthermore, through the use of ML algorithms, linear discriminant analysis, and random forest, we achieved remarkable accuracy, as high as 0.99, in discerning between strong (cognate) and weak (noncognate) binders. The presented ML approach encompasses easily transferable input features, enabling its broad application to any interactome while facilitating the identification of pivotal residues critical for binding interactions. The predictive power of the generated model was probed on similar protein families from 13 diverse species. Our ML model exhibited commendable performance on these additional data sets, showcasing its reliability and robustness across the species barrier.

Computational Investigation of the Covalent Inhibition Mechanism of Bruton's Tyrosine Kinase by Ibrutinib.

Covalent inhibitors represent a promising class of therapeutic compounds. Nonetheless, rationally designing covalent inhibitors to achieve a right balance between selectivity and reactivity remains extremely challenging. To better understand the covalent binding mechanism, a computational study is carried out using the irreversible covalent inhibitor of Bruton tyrosine kinase (BTK) ibrutinib as an example. A multi-µs classical molecular dynamics trajectory of the unlinked inhibitor is generated to explore the fluctuations of the compound associated with the kinase binding pocket. Then, the reaction pathway leading to the formation of the covalent bond with the cysteine residue at position 481 via a Michael addition is determined using the string method in collective variables on the basis of hybrid quantum mechanical-molecular mechanical (QM/MM) simulations. The reaction pathway shows a strong correlation between the covalent bond formation and the protonation/deprotonation events taking place sequentially in the covalent inhibition reaction, consistent with a 3-step reaction with transient thiolate and enolates intermediate states. Two possible atomistic mechanisms affecting deprotonation/protonation events from the thiolate to the enolate intermediate were observed: a highly correlated direct pathway involving proton transfer to the Cα of the acrylamide warhead from the cysteine involving one or a few water molecules and a more indirect pathway involving a long-lived enolate intermediate state following the escape of the proton to the bulk solution. The results are compared with experiments by simulating the long-time kinetics of the reaction using kinetic modeling.

Folding and misfolding of potassium channel monomers during assembly and tetramerization.

The dynamics and folding of potassium channel pore domain monomers are connected to the kinetics of tetramer assembly. In all-atom molecular dynamics simulations of Kv1.2 and KcsA channels, monomers adopt multiple nonnative conformations while the three helices remain folded. Consistent with this picture, NMR studies also find the monomers to be dynamic and structurally heterogeneous. However, a KcsA construct with a disulfide bridge engineered between the two transmembrane helices has an NMR spectrum with well-dispersed peaks, suggesting that the monomer can be locked into a native-like conformation that is similar to that observed in the folded tetramer. During tetramerization, fluoresence resonance energy transfer (FRET) data indicate that monomers rapidly oligomerize upon insertion into liposomes, likely forming a protein-dense region. Folding within this region occurs along separate fast and slow routes, with τfold ∼40 and 1,500 s, respectively. In contrast, constructs bearing the disulfide bond mainly fold via the faster pathway, suggesting that maintaining the transmembrane helices in their native orientation reduces misfolding. Interestingly, folding is concentration independent despite the tetrameric nature of the channel, indicating that the rate-limiting step is unimolecular and occurs after monomer association in the protein-dense region. We propose that the rapid formation of protein-dense regions may help with the assembly of multimeric membrane proteins by bringing together the nascent components prior to assembly. Finally, despite its name, the addition of KcsA's C-terminal "tetramerization" domain does not hasten the kinetics of tetramerization.

Transition rate theory, spectral analysis, and reactive paths.

The kinetics of a dynamical system dominated by two metastable states is examined from the perspective of the activated-dynamics reactive flux formalism, Markov state eigenvalue spectral decomposition, and committor-based transition path theory. Analysis shows that the different theoretical formulations are consistent, clarifying the significance of the inherent microscopic lag-times that are implicated, and that the most meaningful one-dimensional reaction coordinate in the region of the transition state is along the gradient of the committor in the multidimensional subspace of collective variables. It is shown that the familiar reactive flux activated dynamics formalism provides an effective route to calculate the transition rate in the case of a narrow sharp barrier but much less so in the case of a broad flat barrier. In this case, the standard reactive flux correlation function decays very slowly to the plateau value that corresponds to the transmission coefficient. Treating the committor function as a reaction coordinate does not alleviate all issues caused by the slow relaxation of the reactive flux correlation function. A more efficient activated dynamics simulation algorithm may be achieved from a modified reactive flux weighted by the committor. Simulation results on simple systems are used to illustrate the various conceptual points.

Synthesis, Characterization, and Simulation of Four-Armed Megamolecules.

This paper describes the synthesis, characterization, and modeling of a series of molecules having four protein domains attached to a central core. The molecules were assembled with the "megamolecule" strategy, wherein enzymes react with their covalent inhibitors that are substituted on a linker. Three linkers were synthesized, where each had four oligo(ethylene glycol)-based arms terminated in a para-nitrophenyl phosphonate group that is a covalent inhibitor for cutinase. This enzyme is a serine hydrolase and reacts efficiently with the phosphonate to give a new ester linkage at the Ser-120 residue in the active site of the enzyme. Negative-stain transmission electron microscopy (TEM) images confirmed the architecture of the four-armed megamolecules. These cutinase tetramers were also characterized by X-ray crystallography, which confirmed the active-site serine-phosphonate linkage by electron-density maps. Molecular dynamics simulations of the tetracutinase megamolecules using three different force field setups were performed and compared with the TEM observations. Using the Amberff99SB-disp + pH7 force field, the two-dimensional projection distances of the megamolecules were found to agree with the measured dimensions from TEM. The study described here, which combines high-resolution characterization with molecular dynamics simulations, will lead to a comprehensive understanding of the molecular structures and dynamics for this new class of molecules.

Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields.

Many applications in chemistry, biology, and energy storage/conversion research rely on molecular simulations to provide fundamental insight into structural and transport properties of materials with high ionic concentrations. Whether the system is comprised entirely of ions, like ionic liquids, or is a mixture of a polar solvent with a salt, e.g., liquid electrolytes for battery applications, the presence of ions in these materials results in strong local electric fields polarizing solvent molecules and large ions. To predict properties of such systems from molecular simulations often requires either explicit or mean-field inclusion of the influence of polarization on electrostatic interactions. In this manuscript, we review the pros and cons of different treatments of polarization ranging from the mean-field approaches to the most popular explicit polarization models in molecular dynamics simulations of ionic materials. For each method, we discuss their advantages and disadvantages and emphasize key assumptions as well as their adjustable parameters. Strategies for the development of polarizable models are presented with a specific focus on extracting atomic polarizabilities. Finally, we compare simulations using polarizable and nonpolarizable models for several classes of ionic systems, discussing the underlying physics that each approach includes or ignores, implications for implementation and computational efficiency, and the accuracy of properties predicted by these methods compared to experiments.

String Method with Swarms-of-Trajectories, Mean Drifts, Lag Time, and Committor.

The kinetics of a dynamical system comprising two metastable states is formulated in terms of a finite-time propagator in phase space (position and velocity) adapted to the underdamped Langevin equation. Dimensionality reduction to a subspace of collective variables yields familiar expressions for the propagator, committor, and steady-state flux. A quadratic expression for the steady-state flux between the two metastable states can serve as a robust variational principle to determine an optimal approximate committor expressed in terms of a set of collective variables. The theoretical formulation is exploited to clarify the foundation of the string method with swarms-of-trajectories, which relies on the mean drift of short trajectories to determine the optimal transition pathway. It is argued that the conditions for Markovity within a subspace of collective variables may not be satisfied with an arbitrary short time-step and that proper kinetic behaviors appear only when considering the effective propagator for longer lag times. The effective propagator with finite lag time is amenable to an eigenvalue-eigenvector spectral analysis, as elaborated previously in the context of position-based Markov models. The time-correlation functions calculated by swarms-of-trajectories along the string pathway constitutes a natural extension of these developments. The present formulation provides a powerful theoretical framework to characterize the optimal pathway between two metastable states of a system.

A critical perspective on Markov state model treatments of protein-protein association using coarse-grained simulations.

Atomic-level information is essential to explain the specific interactions governing protein-protein recognition in terms of structure and dynamics. Of particular interest is a characterization of the time-dependent kinetic aspects of protein-protein association and dissociation. A powerful framework to characterize the dynamics of complex molecular systems is provided by Markov State Models (MSMs). The central idea is to construct a reduced stochastic model of the full system by defining a set of conformational featured microstates and determining the matrix of transition probabilities between them. While a MSM framework can sometimes be very effective, different combinations of input featurization and simulation methods can significantly affect the robustness and the quality of the information generated from MSMs in the context of protein association. Here, a systematic examination of a variety of MSMs methodologies is undertaken to clarify these issues. To circumvent the uncertainties caused by sampling issues, we use a simplified coarse-grained model of the barnase-barstar protein complex. A sensitivity analysis is proposed to identify the microstates of an MSM that contribute most to the error in conjunction with the transition-based reweighting analysis method for a more efficient and accurate MSM construction.

Engineering Dirhodium Artificial Metalloenzymes for Diazo Coupling Cascade Reactions*.

Artificial metalloenzymes (ArMs) are commonly used to control the stereoselectivity of catalytic reactions, but controlling chemoselectivity remains challenging. In this study, we engineer a dirhodium ArM to catalyze diazo cross-coupling to form an alkene that, in a one-pot cascade reaction, is reduced to an alkane with high enantioselectivity (typically >99 % ee) by an alkene reductase. The numerous protein and small molecule components required for the cascade reaction had minimal effect on ArM catalysis. Directed evolution of the ArM led to improved yields and E/Z selectivities for a variety of substrates, which translated to cascade reaction yields. MD simulations of ArM variants were used to understand the structural role of the cofactor on ArM conformational dynamics. These results highlight the ability of ArMs to control both catalyst stereoselectivity and chemoselectivity to enable reactions in complex media that would otherwise lead to undesired side reactions.

CHARMM-GUI DEER facilitator for spin-pair distance distribution calculations and preparation of restrained-ensemble molecular dynamics simulations.

The double electron-electron resonance (DEER) is a powerful structural biology technique to obtain distance information in the range of 18 to 80 å by measuring the dipolar coupling between two unpaired electron spins. The distance distributions obtained from the experiment provide valuable structural information about the protein in its native environment that can be exploited using restrained ensemble molecular dynamics (reMD) simulations. We present a new tool DEER Facilitator in CHARMM-GUI that consists of two modules Spin-Pair Distributor and reMD Prepper to setup simulations that utilize information from DEER experiments. Spin-Pair Distributor provides a web-based interface to calculate the spin-pair distance distribution of labeled sites in a protein using MD simulations. The calculated distribution can be used to guide the selection of the labeling sites in experiments as well as validate different protein structure models. reMD Prepper facilities the setup of reMD simulations using different types of spin labels in four different environments including vacuum, solution, micelle, and bilayer. The applications of these two modules are demonstrated with several test cases. Spin-Pair Distributor and reMD Prepper are available at http://www.charmm-gui.org/input/deer and http://www.charmm-gui.org/input/deerre. DEER Facilitator is expected to facilitate advanced biomolecular modeling and simulation, thereby leading to an improved understanding of the structure and dynamics of complex biomolecular systems based on experimental DEER data. © 2019 Wiley Periodicals, Inc.

Quantum Chemical Methods for Modeling Covalent Modification of Biological Thiols.

Targeted covalent inhibitor drugs require computational methods that go beyond simple molecular-mechanical force fields in order to model the chemical reactions that occur when they bind to their targets. Here, several semiempirical and density-functional theory (DFT) methods are assessed for their ability to describe the potential energy surface and reaction energies of the covalent modification of a thiol by an electrophile. Functionals such as PBE and B3LYP fail to predict a stable enolate intermediate. This is largely due to delocalization error, which spuriously stabilizes the prereaction complex, in which excess electron density is transferred from the thiolate to the electrophile. Functionals with a high-exact exchange component, range-separated DFT functionals, and variationally optimized exact exchange (i.e., the LC-B05minV functional) correct this issue to various degrees. The large gradient behavior of the exchange enhancement factor is also found to significantly affect the results, leading to the improved performance of PBE0. While ωB97X-D and M06-2X were reasonably accurate, no method provided quantitative accuracy for all three electrophiles, making this a very strenuous test of functional performance. Additionally, one drawback of M06-2X was that molecular dynamics (MD) simulations using this functional were only stable if a fine integration grid was used. The low-cost semiempirical methods, PM3, AM1, and PM7, provide a qualitatively correct description of the reaction mechanism, although the energetics is not quantitatively reliable. As a proof of concept, the potential of mean force for the addition of methylthiolate to methylvinyl ketone was calculated using quantum mechanical/molecular mechanical MD in an explicit polarizable aqueous solvent. © 2019 Wiley Periodicals, Inc.

Statistical mechanics of polarizable force fields based on classical Drude oscillators with dynamical propagation by the dual-thermostat extended Lagrangian.

Polarizable force fields based on classical Drude oscillators offer a practical and computationally efficient avenue to carry out molecular dynamics (MD) simulations of large biomolecular systems. To treat the polarizable electronic degrees of freedom, the Drude model introduces a virtual charged particle that is attached to its parent nucleus via a harmonic spring. Traditionally, the need to relax the electronic degrees of freedom for each fixed set of nuclear coordinates is achieved by performing an iterative self-consistent field (SCF) calculation to satisfy a selected tolerance. This is a computationally demanding procedure that can increase the computational cost of MD simulations by nearly one order of magnitude. To avoid the costly SCF procedure, a small mass is assigned to the Drude particles, which are then propagated as dynamic variables during the simulations via a dual-thermostat extended Lagrangian algorithm. To help clarify the significance of the dual-thermostat extended Lagrangian propagation in the context of the polarizable force field based on classical Drude oscillators, the statistical mechanics of a dual-temperature canonical ensemble is formulated. The conditions for dynamically maintaining the dual-temperature properties in the case of the classical Drude oscillator are analyzed using the generalized Langevin equation.

Scalable molecular dynamics on CPU and GPU architectures with NAMD.

NAMDis a molecular dynamics program designed for high-performance simulations of very large biological objects on CPU- and GPU-based architectures. NAMD offers scalable performance on petascale parallel supercomputers consisting of hundreds of thousands of cores, as well as on inexpensive commodity clusters commonly found in academic environments. It is written in C++ and leans on Charm++ parallel objects for optimal performance on low-latency architectures. NAMD is a versatile, multipurpose code that gathers state-of-the-art algorithms to carry out simulations in apt thermodynamic ensembles, using the widely popular CHARMM, AMBER, OPLS, and GROMOS biomolecular force fields. Here, we review the main features of NAMD that allow both equilibrium and enhanced-sampling molecular dynamics simulations with numerical efficiency. We describe the underlying concepts utilized by NAMD and their implementation, most notably for handling long-range electrostatics; controlling the temperature, pressure, and pH; applying external potentials on tailored grids; leveraging massively parallel resources in multiple-copy simulations; and hybrid quantum-mechanical/molecular-mechanical descriptions. We detail the variety of options offered by NAMD for enhanced-sampling simulations aimed at determining free-energy differences of either alchemical or geometrical transformations and outline their applicability to specific problems. Last, we discuss the roadmap for the development of NAMD and our current efforts toward achieving optimal performance on GPU-based architectures, for pushing back the limitations that have prevented biologically realistic billion-atom objects to be fruitfully simulated, and for making large-scale simulations less expensive and easier to set up, run, and analyze. NAMD is distributed free of charge with its source code at www.ks.uiuc.edu.

Chemical substitutions in the selectivity filter of potassium channels do not rule out constricted-like conformations for C-type inactivation.

In many K+ channels, prolonged activating stimuli lead to a time-dependent reduction in ion conduction, a phenomenon known as C-type inactivation. X-ray structures of the KcsA channel suggest that this inactivated state corresponds to a "constricted" conformation of the selectivity filter. However, the functional significance of the constricted conformation has become a matter of debate. Functional and structural studies based on chemically modified semisynthetic KcsA channels along the selectivity filter led to the conclusion that the constricted conformation does not correspond to the C-type inactivated state. The main results supporting this view include the observation that C-type inactivation is not suppressed by a substitution of D-alanine at Gly77, even though this modification is believed to lock the selectivity filter into its conductive conformation, whereas it is suppressed following amide-to-ester backbone substitutions at Gly77 and Tyr78, even though these structure-conserving modifications are not believed to prevent the selectivity filter from adopting the constricted conformation. However, several untested assumptions about the structural and functional impact of these chemical modifications underlie these arguments. To make progress, molecular dynamics simulations based on atomic models of the KcsA channel were performed. The computational results support the notion that the constricted conformation of the selectivity filter corresponds to the functional C-type inactivated state of the KcsA. Importantly, MD simulations reveal that the semisynthetic KcsAD-ala77 channel can adopt an asymmetrical constricted-like nonconductive conformation and that the amide-to-ester backbone substitutions at Gly77 and Tyr78 perturb the hydrogen bonding involving the buried water molecules stabilizing the constricted conformation.

Perplexing cooperative folding and stability of a low-sequence complexity, polyproline 2 protein lacking a hydrophobic core.

The burial of hydrophobic side chains in a protein core generally is thought to be the major ingredient for stable, cooperative folding. Here, we show that, for the snow flea antifreeze protein (sfAFP), stability and cooperativity can occur without a hydrophobic core, and without α-helices or ß-sheets. sfAFP has low sequence complexity with 46% glycine and an interior filled only with backbone H-bonds between six polyproline 2 (PP2) helices. However, the protein folds in a kinetically two-state manner and is moderately stable at room temperature. We believe that a major part of the stability arises from the unusual match between residue-level PP2 dihedral angle bias in the unfolded state and PP2 helical structure in the native state. Additional stabilizing factors that compensate for the dearth of hydrophobic burial include shorter and stronger H-bonds, and increased entropy in the folded state. These results extend our understanding of the origins of cooperativity and stability in protein folding, including the balance between solvent and polypeptide chain entropies.

Phosphoantigen-induced conformational change of butyrophilin 3A1 (BTN3A1) and its implication on Vγ9Vδ2 T cell activation.

Human Vγ9Vδ2 T cells respond to microbial infections as well as certain types of tumors. The key initiators of Vγ9Vδ2 activation are small, pyrophosphate-containing molecules called phosphoantigens (pAgs) that are present in infected cells or accumulate intracellularly in certain tumor cells. Recent studies demonstrate that initiation of the Vγ9Vδ2 T cell response begins with sensing of pAg via the intracellular domain of the butyrophilin 3A1 (BTN3A1) molecule. However, it is unknown how downstream events can ultimately lead to T cell activation. Here, using NMR spectrometry and molecular dynamics (MD) simulations, we characterize a global conformational change in the B30.2 intracellular domain of BTN3A1 induced by pAg binding. We also reveal by crystallography two distinct dimer interfaces in the BTN3A1 full-length intracellular domain, which are stable in MD simulations. These interfaces lie in close proximity to the pAg-binding pocket and contain clusters of residues that experience major changes of chemical environment upon pAg binding. This suggests that pAg binding disrupts a preexisting conformation of the BTN3A1 intracellular domain. Using a combination of biochemical, structural, and cellular approaches we demonstrate that the extracellular domains of BTN3A1 adopt a V-shaped conformation at rest, and that locking them in this resting conformation without perturbing their membrane reorganization properties diminishes pAg-induced T cell activation. Based on these results, we propose a model in which a conformational change in BTN3A1 is a key event of pAg sensing that ultimately leads to T cell activation.

Crystal Structure and Conformational Dynamics of Pyrococcus furiosus Prolyl Oligopeptidase.

Enzymes in the prolyl oligopeptidase family possess unique structures and substrate specificities that are important for their biological activity and for potential biocatalytic applications. The crystal structures of Pyrococcus furiosus ( Pfu) prolyl oligopeptidase (POP) and the corresponding S477C mutant were determined to 1.9 and 2.2 Å resolution, respectively. The wild type enzyme crystallized in an open conformation, indicating that this state is readily accessible, and it contained bound chloride ions and a prolylproline ligand. These structures were used as starting points for molecular dynamics simulations of Pfu POP conformational dynamics. The simulations showed that large-scale domain opening and closing occurred spontaneously, providing facile substrate access to the active site. Movement of the loop containing the catalytically essential histidine into a conformation similar to those found in structures with fully formed catalytic triads also occurred. This movement was modulated by chloride binding, providing a rationale for experimentally observed activation of POP peptidase catalysis by chloride. Thus, the structures and simulations reported in this study, combined with existing biochemical data, provide a number of insights into POP catalysis.