"Martinizing" the Variational Implicit Solvent Method (VISM): Solvation Free Energy for Coarse-Grained Proteins

Solvation is a fundamental driving force in many biological processes including biomolecular recognition and self-assembly, not to mention protein folding, dynamics, and function. The variational implicit solvent method (VISM) is a theoretical tool currently developed and optimized to estimate solvation free energies for systems of very complex topology, such as biomolecules. VISM's theoretical framework makes it unique because it couples hydrophobic, van der Waals, and electrostatic interactions as a functional of the solvation interface. By minimizing this functional, VISM produces the solvation interface as an output of the theory. In this work, we push VISM to larger scale applications by combining it with coarse-grained solute Hamiltonians adapted from the MARTINI framework, a well-established mesoscale force field for modeling large-scale biomolecule assemblies. We show how MARTINI-VISM (^MVISM) compares with atomistic VISM (^AVISM) for a small set of proteins differing in size, shape, and charge distribution. We also demonstrate ^MVISM's suitability to study the solvation properties of an interesting encounter complex, barnase-barstar. The promising results suggest that coarse-graining the protein with the MARTINI force field is indeed a valuable step to broaden VISM's and MARTINI's applications in the near future.

Spectroscopic and Computational Investigations of Ligand Binding to IspH: Discovery of Non-diphosphate Inhibitors

Isoprenoid biosynthesis is an important area for anti-infective drug development. One isoprenoid target is (E)-1-hydroxy-2-methyl-but-2-enyl 4-diphosphate (HMBPP) reductase (IspH), which forms isopentenyl diphosphate and dimethylallyl diphosphate from HMBPP in a 2H+ /2e- reduction. IspH contains a 4 Fe-4 S cluster, and in this work, we first investigated how small molecules bound to the cluster by using HYSCORE and NRVS spectroscopies. The results of these, as well as other structural and spectroscopic investigations, led to the conclusion that, in most cases, ligands bound to IspH 4 Fe-4 S clusters by η1 coordination, forming tetrahedral geometries at the unique fourth Fe, ligand side chains preventing further ligand (e.g., H2 O, O2 ) binding. Based on these ideas, we used in silico methods to find drug-like inhibitors that might occupy the HMBPP substrate binding pocket and bind to Fe, leading to the discovery of a barbituric acid analogue with a Ki value of ≈500 nm against Pseudomonas aeruginosa IspH.

Activation mechanisms of the first sphingosine-1-phosphate receptor

Activation of the first sphingosine-1-phosphate receptor (S1PR₁ ) promotes permeability of the blood brain barrier, astrocyte and neuronal protection, and lymphocyte egress from secondary lymphoid tissues. Although an agonist often activates the S1PR₁ , the receptor exhibits high levels of basal activity. In this study, we performed long-timescale molecular dynamics and accelerated molecular dynamics (aMD) simulations to investigate activation mechanisms of the ligand-free (apo) S1PR₁ . In the aMD enhanced sampling simulations, we observed four independent events of activation, which is characterized by close interaction between Y311^7.53 and Y221^5.58 and increased distance between the intracellular ends of transmembrane (TM) helices 3 and 6. Although TM helices TM3, TM6, TM5 and, TM7 are associated with GPCR activation, we discovered that their movements are not necessarily correlated during activation. Instead, TM5 showed a decreased correlation with each of these regions during activation. During activation of the apo receptor, Y221^5.58 and Y311^7.53 became more solvated, because a water channel formed in the intracellular pocket. Additionally, a lipid molecule repeatedly entered the receptor between the extracellular ends of TM1 and TM7, providing important insights into the pathway of ligand entry into the S1PR₁ .

Effect of donor atom identity on metal-binding pharmacophore coordination

The inhibition and binding of three metal-binding pharmacophores (MBPs), 2-hydroxycyclohepta-2,4,6-trien-1-one (tropolone), 2-mercaptopyridine-N-oxide (1,2-HOPTO), and 2-hydroxycyclohepta-2,4,6-triene-1-thione (thiotropolone) to human carbonic anhydrase II (hCAII) and a mutant protein hCAII L198G were investigated. These MBPs displayed bidentate coordination to the active site Zn(II) metal ion, but the MBPs respond to the mutation of L198G differently, as characterized by inhibition activity assays and X-ray crystallography. The L198G mutation increases the active site volume thereby decreasing the steric pressure exerted on MBPs upon binding, allowing changes in MBP coordination to be observed. When comparing the binding mode of tropolone to thiotropolone or 1,2-HOPTO (O,O versus O,S donor sets), structural modifications of the hCAII active site were shown to have a stronger effect on MBPs with an O,O versus O,S donor set. These findings were corroborated with density functional theory (DFT) calculations of model coordination complexes. These results suggest that the MBP binding geometry is a malleable interaction, particularly for certain ligands, and that the identity of the donor atoms influences the response of the ligand to changes in the protein active site environment. Understanding underlying interactions between a MBP and a metalloenzyme active site may aid in the design and development of potent metalloenzyme inhibitors.

NMR structure-based optimization of Staphylococcus aureus sortase A pyridazinone inhibitors

Staphylococcus aureus is a leading cause of hospital-acquired infections in the USA and is a major health concern as methicillin-resistant S. aureus and other antibiotic-resistant strains are common. Compounds that inhibit the S. aureus sortase (SrtA) cysteine transpeptidase may function as potent anti-infective agents as this enzyme attaches virulence factors to the bacterial cell wall. While a variety of SrtA inhibitors have been discovered, the vast majority of these small molecules have not been optimized using structure-based approaches. Here we have used NMR spectroscopy to determine the molecular basis through which pyridazinone-based small molecules inhibit SrtA. These inhibitors covalently modify the active cysteine thiol and partially mimic the natural substrate of SrtA by inducing the closure of an active site loop. Computational and synthetic chemistry methods led to second-generation analogues that are ~70-fold more potent than the lead molecule. These optimized molecules exhibit broad-spectrum activity against other types of class A sortases, have reduced cytotoxicity, and impair SrtA-mediated protein display on S. aureus cell surface. Our work shows that pyridazinone analogues are attractive candidates for further development into anti-infective agents, and highlights the utility of employing NMR spectroscopy and solubility-optimized small molecules in structure-based drug discovery.

Gaussian Accelerated Molecular Dynamics in NAMD

Gaussian accelerated molecular dynamics (GaMD) is a recently developed enhanced sampling technique that provides efficient free energy calculations of biomolecules. Like the previous accelerated molecular dynamics (aMD), GaMD allows for “unconstrained” enhanced sampling without the need to set predefined collective variables and so is useful for studying complex biomolecular conformational changes such as protein folding and ligand binding. Furthermore, because the boost potential is constructed using a harmonic function that follows Gaussian distribution in GaMD, cumulant expansion to the second order can be applied to recover the original free energy profiles of proteins and other large biomolecules, which solves a long-standing energetic reweighting problem of the previous aMD method. Taken together, GaMD offers major advantages for both unconstrained enhanced sampling and free energy calculations of large biomolecules. Here, we have implemented GaMD in the NAMD package on top of the existing aMD feature and validated it on three model systems: alanine dipeptide, the chignolin fast-folding protein, and the M₃ muscarinic G protein-coupled receptor (GPCR). For alanine dipeptide, while conventional molecular dynamics (cMD) simulations performed for 30 ns are poorly converged, GaMD simulations of the same length yield free energy profiles that agree quantitatively with those of 1000 ns cMD simulation. Further GaMD simulations have captured folding of the chignolin and binding of the acetylcholine (ACh) endogenous agonist to the M₃ muscarinic receptor. The reweighted free energy profiles are used to characterize the protein folding and ligand binding pathways quantitatively. GaMD implemented in the scalable NAMD is widely applicable to enhanced sampling and free energy calculations of large biomolecules.

Striking Plasticity of CRISPR-Cas9 and Key Role of Non-target DNA, as Revealed by Molecular Simulations

The CRISPR (clustered regularly interspaced short palindromic repeats)-Cas9 system recently emerged as a transformative genome-editing technology that is innovating basic bioscience and applied medicine and biotechnology. The endonuclease Cas9 associates with a guide RNA to match and cleave complementary sequences in double stranded DNA, forming an RNA:DNA hybrid and a displaced non-target DNA strand. Although extensive structural studies are ongoing, the conformational dynamics of Cas9 and its interplay with the nucleic acids during association and DNA cleavage are largely unclear. Here, by employing multi-microsecond time scale molecular dynamics, we reveal the conformational plasticity of Cas9 and identify key determinants that allow its large-scale conformational changes during nucleic acid binding and processing. We show how the "closure" of the protein, which accompanies nucleic acid binding, fundamentally relies on highly coupled and specific motions of the protein domains, collectively initiating the prominent conformational changes needed for nucleic acid association. We further reveal a key role of the non-target DNA during the process of activation of the nuclease HNH domain, showing how the nontarget DNA positioning triggers local conformational changes that favor the formation of a catalytically competent Cas9. Finally, a remarkable conformational plasticity is identified as an intrinsic property of the HNH domain, constituting a necessary element that allows for the HNH repositioning. These novel findings constitute a reference for future experimental studies aimed at a full characterization of the dynamic features of the CRISPR-Cas9 system, and-more importantly-call for novel structure engineering efforts that are of fundamental importance for the rational design of new genome-engineering applications.

Dynamic Structure and Inhibition of a Malaria Drug Target: Geranylgeranyl Diphosphate Synthase

We report a molecular dynamics investigation of the structure, function, and inhibition of geranylgeranyl diphosphate synthase (GGPPS), a potential drug target, from the malaria parasite Plasmodium vivax. We discovered several GGPPS inhibitors, benzoic acids, and determined their structures crystallographically. We then used molecular dynamics simulations to investigate the dynamics of three such inhibitors and two bisphosphonate inhibitors, zoledronate and a lipophilic analogue of zoledronate, as well as the enzyme's product, GGPP. We were able to identify the main motions that govern substrate binding and product release as well as the molecular features required for GGPPS inhibition by both classes of inhibitor. The results are of broad general interest because they represent the first detailed investigation of the mechanism of action, and inhibition, of an important antimalarial drug target, geranylgeranyl diphosphate synthase, and may help guide the development of other, novel inhibitors as new drug leads.

Stochastic level-set variational implicit-solvent approach to solute-solvent interfacial fluctuations

Recent years have seen the initial success of a variational implicit-solvent model (VISM), implemented with a robust level-set method, in capturing efficiently different hydration states and providing quantitatively good estimation of solvation free energies of biomolecules. The level-set minimization of the VISM solvation free-energy functional of all possible solute-solvent interfaces or dielectric boundaries predicts an equilibrium biomolecular conformation that is often close to an initial guess. In this work, we develop a theory in the form of Langevin geometrical flow to incorporate solute-solvent interfacial fluctuations into the VISM. Such fluctuations are crucial to biomolecular conformational changes and binding process. We also develop a stochastic level-set method to numerically implement such a theory. We describe the interfacial fluctuation through the "normal velocity" that is the solute-solvent interfacial force, derive the corresponding stochastic level-set equation in the sense of Stratonovich so that the surface representation is independent of the choice of implicit function, and develop numerical techniques for solving such an equation and processing the numerical data. We apply our computational method to study the dewetting transition in the system of two hydrophobic plates and a hydrophobic cavity of a synthetic host molecule cucurbit[7]uril. Numerical simulations demonstrate that our approach can describe an underlying system jumping out of a local minimum of the free-energy functional and can capture dewetting transitions of hydrophobic systems. In the case of two hydrophobic plates, we find that the wavelength of interfacial fluctuations has a strong influence to the dewetting transition. In addition, we find that the estimated energy barrier of the dewetting transition scales quadratically with the inter-plate distance, agreeing well with existing studies of molecular dynamics simulations. Our work is a first step toward the inclusion of fluctuations into the VISM and understanding the impact of interfacial fluctuations on biomolecular solvation with an implicit-solvent approach.

Unconstrained Enhanced Sampling for Free Energy Calculations of Biomolecules: A Review

Free energy calculations are central to understanding the structure, dynamics and function of biomolecules. Yet insufficient sampling of biomolecular configurations is often regarded as one of the main sources of error. Many enhanced sampling techniques have been developed to address this issue. Notably, enhanced sampling methods based on biasing collective variables (CVs), including the widely used umbrella sampling, adaptive biasing force and metadynamics, have been discussed in a recent excellent review (Abrams and Bussi, Entropy, 2014). Here, we aim to review enhanced sampling methods that do not require predefined system-dependent CVs for biomolecular simulations and as such do not suffer from the hidden energy barrier problem as encountered in the CV-biasing methods. These methods include, but are not limited to, replica exchange/parallel tempering, self-guided molecular/Langevin dynamics, essential energy space random walk and accelerated molecular dynamics. While it is overwhelming to describe all details of each method, we provide a summary of the methods along with the applications and offer our perspectives. We conclude with challenges and prospects of the unconstrained enhanced sampling methods for accurate biomolecular free energy calculations.

Improved cryoEM-Guided Iterative Molecular Dynamics--Rosetta Protein Structure Refinement Protocol for High Precision Protein Structure Prediction

Many excellent methods exist that incorporate cryo-electron microscopy (cryoEM) data to constrain computational protein structure prediction and refinement. Previously, it was shown that iteration of two such orthogonal sampling and scoring methods – Rosetta and molecular dynamics (MD) simulations – facilitated exploration of conformational space in principle. Here, we go beyond a proof-of-concept study and address significant remaining limitations of the iterative MD–Rosetta protein structure refinement protocol. Specifically, all parts of the iterative refinement protocol are now guided by medium-resolution cryoEM density maps, and previous knowledge about the native structure of the protein is no longer necessary. Models are identified solely based on score or simulation time. All four benchmark proteins showed substantial improvement through three rounds of the iterative refinement protocol. The best-scoring final models of two proteins had sub-Ångstrom RMSD to the native structure over residues in secondary structure elements. Molecular dynamics was most efficient in refining secondary structure elements and was thus highly complementary to the Rosetta refinement which is most powerful in refining side chains and loop regions.

How to deal with multiple binding poses in alchemical relative protein-ligand binding free energy calculations

Recent advances in improved force fields and sampling methods have made it possible for the accurate calculation of protein–ligand binding free energies. Alchemical free energy perturbation (FEP) using an explicit solvent model is one of the most rigorous methods to calculate relative binding free energies. However, for cases where there are high energy barriers separating the relevant conformations that are important for ligand binding, the calculated free energy may depend on the initial conformation used in the simulation due to the lack of complete sampling of all the important regions in phase space. This is particularly true for ligands with multiple possible binding modes separated by high energy barriers, making it difficult to sample all relevant binding modes even with modern enhanced sampling methods. In this paper, we apply a previously developed method that provides a corrected binding free energy for ligands with multiple binding modes by combining the free energy results from multiple alchemical FEP calculations starting from all enumerated poses, and the results are compared with Glide docking and MM-GBSA calculations. From these calculations, the dominant ligand binding mode can also be predicted. We apply this method to a series of ligands that bind to c-Jun N-terminal kinase-1 (JNK1) and obtain improved free energy results. The dominant ligand binding modes predicted by this method agree with the available crystallography, while both Glide docking and MM-GBSA calculations incorrectly predict the binding modes for some ligands. The method also helps separate the force field error from the ligand sampling error, such that deviations in the predicted binding free energy from the experimental values likely indicate possible inaccuracies in the force field. An error in the force field for a subset of the ligands studied was identified using this method, and improved free energy results were obtained by correcting the partial charges assigned to the ligands. This improved the root-mean-square error (RMSE) for the predicted binding free energy from 1.9 kcal/mol with the original partial charges to 1.3 kcal/mol with the corrected partial charges.

Identification of protein-ligand binding sites by the level-set variational implicit-solvent approach

Protein–ligand binding is a key biological process at the molecular level. The identification and characterization of small-molecule binding sites on therapeutically relevant proteins have tremendous implications for target evaluation and rational drug design. In this work, we used the recently developed level-set variational implicit-solvent model (VISM) with the Coulomb field approximation (CFA) to locate and characterize potential protein–small-molecule binding sites. We applied our method to a data set of 515 protein–ligand complexes and found that 96.9% of the cocrystallized ligands bind to the VISM-CFA-identified pockets and that 71.8% of the identified pockets are occupied by cocrystallized ligands. For 228 tight-binding protein–ligand complexes (i.e, complexes with experimental pK_d values larger than 6), 99.1% of the cocrystallized ligands are in the VISM-CFA-identified pockets. In addition, it was found that the ligand binding orientations are consistent with the hydrophilic and hydrophobic descriptions provided by VISM. Quantitative characterization of binding pockets with topological and physicochemical parameters was used to assess the “ligandability” of the pockets. The results illustrate the key interactions between ligands and receptors and can be very informative for rational drug design.

Enzyme localization, crowding, and buffers collectively modulate diffusion-influenced signal transduction: Insights from continuum diffusion modeling

Biochemical reaction networks consisting of coupled enzymes connect substrate signaling events with biological function. Substrates involved in these reactions can be strongly influenced by diffusion "barriers" arising from impenetrable cellular structures and macromolecules, as well as interactions with biomolecules, especially within crowded environments. For diffusion-influenced reactions, the spatial organization of diffusion barriers arising from intracellular structures, non-specific crowders, and specific-binders (buffers) strongly controls the temporal and spatial reaction kinetics. In this study, we use two prototypical biochemical reactions, a Goodwin oscillator, and a reaction with a periodic source/sink term to examine how a diffusion barrier that partitions substrates controls reaction behavior. Namely, we examine how conditions representative of a densely packed cytosol, including reduced accessible volume fraction, non-specific interactions, and buffers, impede diffusion over nanometer length-scales. We find that diffusion barriers can modulate the frequencies and amplitudes of coupled diffusion-influenced reaction networks, as well as give rise to "compartments" of decoupled reactant populations. These effects appear to be intensified in the presence of buffers localized to the diffusion barrier. These findings have strong implications for the role of the cellular environment in tuning the dynamics of signaling pathways.

Molecular dynamic study of MlaC protein in Gram-negative bacteria: conformational flexibility, solvent effect and protein-phospholipid binding

The composition of the outer membrane in Gram-negative bacteria is asymmetric, with the lipopolysaccharides found in the outer leaflet and phospholipids in the inner leaflet. The MlaC protein transfers phospholipids from the outer to inner membrane to maintain such lipid asymmetry in the Mla pathway. In this work, we have performed molecular dynamics simulations on apo and phospholipid-bound systems to study the dynamical properties of MlaC. Our simulations show that the phospholipid forms hydrophobic interactions with the protein. Residues surrounding the entrance of the binding site exhibit correlated motions to control the site opening and closing. Lipid binding leads to increase of the binding pocket volume and precludes entry of the water molecules. However, in the absence of the phospholipid, water molecules can freely move in and out of the binding site when the pocket is open. Dehydration occurs when the pocket closes. This study provides dynamic information of the MlaC protein and may facilitate the design of antibiotics against the Mla pathway of Gram-negative bacteria.

Development of Potent and Selective Inhibitors for Group VIA Calcium-Independent Phospholipase A2 Guided by Molecular Dynamics and Structure-Activity Relationships

The development of inhibitors for phospholipase A2 (PLA2) is important in elucidating the enzymes implication in various biological pathways. PLA2 enzymes are an important pharmacological target implicated in various inflammatory diseases. Computational chemistry, organic synthesis, and in vitro assays were employed to develop potent and selective inhibitors for group VIA calcium-independent PLA2. A set of fluoroketone inhibitors was studied for their binding mode with two human cytosolic PLA2 enzymes: group IVA cPLA2 and group VIA iPLA2. New compounds were synthesized and assayed toward three major PLA2s. This study led to the development of four potent and selective thioether fluoroketone inhibitors as well as a thioether keto-1,2,4-oxadiazole inhibitor for GVIA iPLA2, which will serve as lead compounds for future development and studies. The keto-1,2,4-oxadiazole functionality with a thioether is a novel structure, and it will be used as a lead to develop inhibitors with higher potency and selectivity toward GVIA iPLA2.

Hybrid finite element and Brownian dynamics method for charged particles

Diffusion is often the rate-determining step in many biological processes. Currently, the two main computational methods for studying diffusion are stochastic methods, such as Brownian dynamics, and continuum methods, such as the finite element method. A previous study introduced a new hybrid diffusion method that couples the strengths of each of these two methods, but was limited by the lack of interactions among the particles; the force on each particle had to be from an external field. This study further develops the method to allow charged particles. The method is derived for a general multidimensional system and is presented using a basic test case for a one-dimensional linear system with one charged species and a radially symmetric system with three charged species.

General trends of dihedral conformational transitions in a globular protein

Dihedral conformational transitions are analyzed systematically in a model globular protein, cytochrome P450cam, to examine their structural and chemical dependences through combined conventional molecular dynamics (cMD), accelerated molecular dynamics (aMD) and adaptive biasing force (ABF) simulations. The aMD simulations are performed at two acceleration levels, using dihedral and dual boost, respectively. In comparison with cMD, aMD samples protein dihedral transitions approximately two times faster on average using dihedral boost, and ∼ 3.5 times faster using dual boost. In the protein backbone, significantly higher dihedral transition rates are observed in the bend, coil, and turn flexible regions, followed by the β bridge and β sheet, and then the helices. Moreover, protein side chains of greater length exhibit higher transition rates on average in the aMD-enhanced sampling. Side chains of the same length (particularly Nχ = 2) exhibit decreasing transition rates with residues when going from hydrophobic to polar, then charged and aromatic chemical types. The reduction of dihedral transition rates is found to be correlated with increasing energy barriers as identified through ABF free energy calculations. These general trends of dihedral conformational transitions provide important insights into the hierarchical dynamics and complex free energy landscapes of functional proteins.

Computation of pH-dependent binding free energies

Protein-ligand binding accompanies changes in the surrounding electrostatic environments of the two binding partners and may lead to changes in protonation upon binding. In cases where the complex formation results in a net transfer of protons, the binding process is pH-dependent. However, conventional free energy computations or molecular docking protocols typically employ fixed protonation states for the titratable groups in both binding partners set a priori, which are identical for the free and bound states. In this review, we draw attention to these important yet largely ignored binding-induced protonation changes in protein-ligand association by outlining physical origins and prevalence of the protonation changes upon binding. Following a summary of various theoretical methods for pKa prediction, we discuss the theoretical framework to examine the pH dependence of protein-ligand binding processes.

A self-consistent phase-field approach to implicit solvation of charged molecules with Poisson-Boltzmann electrostatics

Dielectric boundary based implicit-solvent models provide efficient descriptions of coarse-grained effects, particularly the electrostatic effect, of aqueous solvent. Recent years have seen the initial success of a new such model, variational implicit-solvent model (VISM) [Dzubiella, Swanson, and McCammon Phys. Rev. Lett. 96, 087802 (2006) and J. Chem. Phys. 124, 084905 (2006)], in capturing multiple dry and wet hydration states, describing the subtle electrostatic effect in hydrophobic interactions, and providing qualitatively good estimates of solvation free energies. Here, we develop a phase-field VISM to the solvation of charged molecules in aqueous solvent to include more flexibility. In this approach, a stable equilibrium molecular system is described by a phase field that takes one constant value in the solute region and a different constant value in the solvent region, and smoothly changes its value on a thin transition layer representing a smeared solute-solvent interface or dielectric boundary. Such a phase field minimizes an effective solvation free-energy functional that consists of the solute-solvent interfacial energy, solute-solvent van der Waals interaction energy, and electrostatic free energy described by the Poisson-Boltzmann theory. We apply our model and methods to the solvation of single ions, two parallel plates, and protein complexes BphC and p53/MDM2 to demonstrate the capability and efficiency of our approach at different levels. With a diffuse dielectric boundary, our new approach can describe the dielectric asymmetry in the solute-solvent interfacial region. Our theory is developed based on rigorous mathematical studies and is also connected to the Lum-Chandler-Weeks theory (1999). We discuss these connections and possible extensions of our theory and methods.