Extending the Applicability of the Nonlinear Poisson−Boltzmann Equation: Multiple Dielectric Constants and Multivalent Ions

Implicit solvent methods for free energy estimation

Solvation is a fundamental contribution in many biological processes and especially in molecular binding. Its estimation can be performed by means of several computational approaches. The aim of this review is to give an overview of existing theories and methods to estimate solvent effects giving a specific focus on the category of implicit solvent models and their use in Molecular Dynamics. In many of these models, the solvent is considered as a continuum homogenous medium, while the solute can be represented at the atomic detail and at different levels of theory. Despite their degree of approximation, implicit methods are still widely employed due to their trade-off between accuracy and efficiency. Their derivation is rooted in the statistical mechanics and integral equations disciplines, some of the related details being provided here. Finally, methods that combine implicit solvent models and molecular dynamics simulation, are briefly described.

A multiscale approach to predicting affinity changes in protein-protein interfaces

Substitution mutations in protein-protein interfaces can have a substantial effect on binding, which has consequences in basic and applied biomedical research. Experimental expression, purification, and affinity determination of protein complexes is an expensive and time-consuming means of evaluating the effect of mutations, making a fast and accurate in silico method highly desirable. When the structure of the wild-type complex is known, it is possible to economically evaluate the effect of point mutations with knowledge based potentials, which do not model backbone flexibility, but these have been validated only for single mutants. Substitution mutations tend to induce local conformational rearrangements only. Accordingly, ZEMu (Zone Equilibration of Mutants) flexibilizes only a small region around the site of mutation, then computes its dynamics under a physics-based force field. We validate with 1254 experimental mutants (with 1-15 simultaneous substitutions) in a wide variety of different protein environments (65 protein complexes), and obtain a significant improvement in the accuracy of predicted ΔΔG.

A combined MPI-CUDA parallel solution of linear and nonlinear Poisson-Boltzmann equation

The Poisson-Boltzmann equation models the electrostatic potential generated by fixed charges on a polarizable solute immersed in an ionic solution. This approach is often used in computational structural biology to estimate the electrostatic energetic component of the assembly of molecular biological systems. In the last decades, the amount of data concerning proteins and other biological macromolecules has remarkably increased. To fruitfully exploit these data, a huge computational power is needed as well as software tools capable of exploiting it. It is therefore necessary to move towards high performance computing and to develop proper parallel implementations of already existing and of novel algorithms. Nowadays, workstations can provide an amazing computational power: up to 10 TFLOPS on a single machine equipped with multiple CPUs and accelerators such as Intel Xeon Phi or GPU devices. The actual obstacle to the full exploitation of modern heterogeneous resources is efficient parallel coding and porting of software on such architectures. In this paper, we propose the implementation of a full Poisson-Boltzmann solver based on a finite-difference scheme using different and combined parallel schemes and in particular a mixed MPI-CUDA implementation. Results show great speedups when using the two schemes, achieving an 18.9x speedup using three GPUs.

Biological applications of classical electrostatics methods

Continuum electrostatics modeling of solvation based on the Poisson–Boltzmann (PB) equation has gained wide acceptance in biomolecular applications such as energetic analysis and structural visualization. Successful application of the PB solvent models requires careful calibration of the solvation parameters. Extensive testing and validation is also important to ensure accuracy in their applications. Limitation in the continuum modeling of solvation is also a known issue in certain biomolecular applications. Growing interest in membrane systems has further spurred developmental efforts to allow inclusion of membrane in the PB solvent models. Despite their past successes due to careful parameterization, algorithm development and parallel implementation, there is still much to be done to improve their transferability from the small molecular systems upon which they were developed and validated to complex macromolecular systems as advances in technology continue to push forward, providing ever greater computational resources to researchers to study more interesting biological systems of higher complexity.

On the Modeling of Polar Component of Solvation Energy using Smooth Gaussian-Based Dielectric Function

Traditional implicit methods for modeling electrostatics in biomolecules use a two-dielectric approach: a biomolecule is assigned low dielectric constant while the water phase is considered as a high dielectric constant medium. However, such an approach treats the biomolecule-water interface as a sharp dielectric border between two homogeneous dielectric media and does not account for inhomogeneous dielectric properties of the macromolecule as well. Recently we reported a new development, a smooth Gaussian-based dielectric function which treats the entire system, the solute and the water phase, as inhomogeneous dielectric medium (J Chem Theory Comput. 2013 Apr 9; 9(4): 2126-2136.). Here we examine various aspects of the modeling of polar solvation energy in such inhomogeneous systems in terms of the solute-water boundary and the inhomogeneity of the solute in the absence of water surrounding. The smooth Gaussian-based dielectric function is implemented in the DelPhi finite-difference program, and therefore the sensitivity of the results with respect to the grid parameters is investigated, and it is shown that the calculated polar solvation energy is almost grid independent. Furthermore, the results are compared with the standard two-media model and it is demonstrated that on average, the standard method overestimates the magnitude of the polar solvation energy by a factor 2.5. Lastly, the possibility of the solute to have local dielectric constant larger than of a bulk water is investigated in a benchmarking test against experimentally determined set of pKa's and it is speculated that side chain rearrangements could result in local dielectric constant larger than 80.

Recent progress in adapting Poisson–Boltzmann methods to molecular simulations

Electrostatic solvation modeling based upon the Poisson–Boltzmann equation is widely used in studies of biomolecular structures and functions. This manuscript provides a thorough review of published efforts to adapt the numerical Poisson–Boltzmann methods to molecular simulations so that these methods can be extended to biomolecular studies involving conformational fluctuation and/or dynamics. We first review the fundamental works on how to define the electrostatic free energy and the Maxwell stress tensor. These topics are followed by three different strategies in developing algorithms to compute electrostatic forces and how to improve their numerical performance. Finally procedures are also presented in detail on how to discretize these algorithms for numerical calculations. Given the pioneer works reviewed here, further developmental efforts will be on how to balance efficiency and accuracy in these theoretical sound approaches — two important issues in applying any numerical algorithms for routine biomolecular applications. Even if not reviewed here, more advanced numerical solvers are certainly necessary to achieve higher accuracy than the widely used classical methods to improve the overall performance of the numerical Poisson–Boltzmann methods.

Accuracy of continuum electrostatic calculations based on three common dielectric boundary definitions

We investigate the influence of three common definitions of the solute/solvent dielectric boundary (DB) on the accuracy of the electrostatic solvation energy ΔG_el computed within the Poisson Boltzmann and the generalized Born models of implicit solvation. The test structures include small molecules, peptides and small proteins; explicit solvent ΔG_el are used as accuracy reference. For common atomic radii sets BONDI, PARSE (and ZAP9 for small molecules) the use of van der Waals (vdW) DB results, on average, in considerably larger errors in ΔG_el than the molecular surface (MS) DB. The optimal probe radius ρ_w for which the MS DB yields the most accurate ΔG_el varies considerably between structure types. The solvent accessible surface (SAS) DB becomes optimal at ρ_w ~ 0.2 Å (exact value is sensitive to the structure and atomic radii), at which point the average accuracy of ΔG_el is comparable to that of the MS-based boundary. The geometric equivalence of SAS to vdW surface based on the same atomic radii uniformly increased by ρ_w gives the corresponding optimal vdW DB. For small molecules, the optimal vdW DB based on BONDI + 0.2 Å radii can yield ΔG_el estimates at least as accurate as those based on the optimal MS DB. Also, in small molecules, pairwise charge-charge interactions computed with the optimal vdW DB are virtually equal to those computed with the MS DB, suggesting that in this case the two boundaries are practically equivalent by the electrostatic energy criteria. In structures other than small molecules, the optimal vdW and MS dielectric boundaries are not equivalent: the respective pairwise electrostatic interactions in the presence of solvent can differ by up to 5 kcal/mol for individual atomic pairs in small proteins, even when the total ΔG_el are equal. For small proteins, the average decrease in pairwise electrostatic interactions resulting from the switch from optimal MS to optimal vdW DB definition can be mimicked within the MS DB definition by doubling of the solute dielectric constant. However, the use of the higher interior dielectric does not eliminate the large individual deviations between pairwise interactions computed within the two DB definitions. It is argued that while the MS based definition of the dielectric boundary is more physically correct in some types of practical calculations, the choice is not so clear in some other common scenarios.

Modified Poisson–Boltzmann equations for characterizing biomolecular solvation

Methods for computing electrostatic interactions often account implicitly for the solvent, due to the much smaller number of degrees of freedom involved. In the Poisson–Boltzmann (PB) approach the electrostatic potential is obtained by solving the Poisson–Boltzmann equation (PBE), where the solvent region is modeled as a homogeneous medium with a high dielectric constant. PB however is not exempt of problems. It does not take into account for example the sizes of the ions in the atmosphere surrounding the solute, nor does it take into account the inhomogeneous dielectric response of water due to the presence of a highly charged surface. In this paper we review two major modifications of PB that circumvent these problems, namely the size-modified PB (SMPB) equation and the Dipolar Poisson–Boltzmann Langevin (DPBL) model. In SMPB, steric effects between ions are accounted for with a lattice gas model. In DPBL, the solvent region is no longer modeled as a homogeneous dielectric media but rather as an assembly of self-orienting interacting dipoles of variable density. This model results in a dielectric profile that transits smoothly from the solute to the solvent region as well as in a variable solvent density that depends on the charges of the solute. We show successful applications of the DPBL formalism to computing the solvation free energies of isolated ions in water. Further developments of more accurately modified PB models are discussed.

Insight into the allosteric inhibition of Abl kinase

Abl kinase inhibitors targeting the ATP binding pocket are currently used as a front-line therapy for the treatment of chronic myelogenous leukemia (CML), but their use has significant limitation because of the development of drug resistance (especially due to the T315I mutation). Two compounds (GNF-2 and BO1) have been found able to inhibit the Abl activity through a peculiar mechanism of action. Particularly, GNF-2 acts as allosteric inhibitor against Bcr-Abl wild type (wt), but it has no activity against the gatekeeper mutant T315I. Its activity against the last mutant reappears when used together with an ATP-competitive inhibitor such as Imatinib or Nilotinib. A crystal structure of GNF-2 bound to the Abl myristoyl pocket (MP) has been released. On the contrary, BO1 shows an ATP-competitive/mixed mechanism of action against the wt, while it acts as an allosteric inhibitor against T315I. In order to better understand the mechanism of Abl allosteric inhibition, MD simulations and MM/GBSA analysis were performed on Abl wt and T315I in complex with GNF-2 and BO1, and the results were compared to those found for the natural myristoyl ligand. Similarly to that observed for the myristoyl group, the binding of an allosteric inhibitor to the MP promotes the formation of a compact and inhibited conformation of the wt protein, characterized by the stabilization of the intramolecular interactions that occur between SH2-SH3 and kinase domains. Conversely, an overall higher flexibility was observed with the Abl T315I mutant, especially in the case of GNF-2. Our analysis highlighted differences in the dynamic behavior of GNF-2 and BO1 which could explain the different biological profiles of the two allosteric inhibitors against the T315I mutant.

Recent host range expansion of canine distemper virus and variation in its receptor, the signaling lymphocyte activation molecule, in carnivores

The signaling lymphocyte activation molecule (SLAM) is a receptor for morbilliviruses. To understand the recent host range expansion of canine distemper virus (CDV) in carnivores, we determined the nucleotide sequences of SLAMs of various carnivores and generated three-dimensional homology SLAM models. Thirty-four amino acid residues were found for the candidates binding to CDV on the interface of the carnivore SLAMs. SLAM of the domestic dog (Canis lupus familiaris) were similar to those of other members of the suborder Caniformia, indicating that the animals in this group have similar sensitivity to dog CDV. However, they were different at nine positions from those of felids. Among the nine residues, four of domestic cat (Felis catus) SLAM (72, 76, 82, and 129) and three of lion (Panthera leo persica) SLAM (72, 82, and 129) were associated with charge alterations, suggesting that the felid interfaces have lower affinities to dog CDV. Only the residue at 76 was different between domestic cat and lion SLAM interfaces. The domestic cat SLAM had threonine at 76, whereas the lion SLAM had arginine, a positively charged residue like that of the dog SLAM. The cat SLAM with threonine is likely to have lower affinity to CDV-H and to confer higher resistance against dog CDV. Thus, the four residues (72, 76, 82, and 129) on carnivore SLAMs are important for the determination of affinity and sensitivity with CDV. Additionally, the CDV-H protein of felid strains had a substitution of histidine for tyrosine at 549 of dog CDV-H and may have higher affinity to lion SLAM. Three-dimensional model construction is a new risk assessment method of morbillivirus infectivity. Because the method is applicable to animals that have no information about virus infection, it is especially useful for morbillivirus risk assessment and wildlife conservation.

Discovery of new sites for drug binding to the hypertension-related renin-angiotensinogen complex

Renin (REN) is a key drug target to stop the hypertension cascade, but thus far only one direct inhibitor has been made commercially available. In this study, we assess an innovative REN inhibition strategy, by targeting the interface of the renin:angiotensinogen (REN:ANG) complex. We characterized the energetic role of interfacial residues of REN:ANG and identified the ones responsible for protein:protein binding, which can serve as drug targets for disruption of the REN:ANG association. For this purpose, we applied a computational alanine scanning mutagenesis protocol, which measures the contribution of each side chain for the protein:protein binding free energy with an accuracy of ≈ 1 kcal/mol. As a result, in REN and ANG, six and eight residues were found to be critical for binding, respectively. The leading force behind REN:ANG complexation was found to be the hydrophobic effect. The binding free energy per residue was found to be proportional to the buried area. Residues responsible for binding were occluded from water at the complex, which promotes an efficient pairing between the two proteins. Two druggable pockets involving critical residues for binding were found on the surface of REN, where small druglike molecules can bind and disrupt the ANG:REN association that may provide an efficient way to achieve REN inhibition and control hypertension.

Electrostatic forces in the Poisson-Boltzmann systems

Continuum modeling of electrostatic interactions based upon numerical solutions of the Poisson-Boltzmann equation has been widely used in structural and functional analyses of biomolecules. A limitation of the numerical strategies is that it is conceptually difficult to incorporate these types of models into molecular mechanics simulations, mainly because of the issue in assigning atomic forces. In this theoretical study, we first derived the Maxwell stress tensor for molecular systems obeying the full nonlinear Poisson-Boltzmann equation. We further derived formulations of analytical electrostatic forces given the Maxwell stress tensor and discussed the relations of the formulations with those published in the literature. We showed that the formulations derived from the Maxwell stress tensor require a weaker condition for its validity, applicable to nonlinear Poisson-Boltzmann systems with a finite number of singularities such as atomic point charges and the existence of discontinuous dielectric as in the widely used classical piece-wise constant dielectric models.

Electrostatics of proteins in dielectric solvent continua. I. An accurate and efficient reaction field description

We present a reaction field (RF) method which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of an electrostatics calculation with polarizable molecular mechanics (MM) force fields. The method combines an approach originally suggested by Egwolf and Tavan [J. Chem. Phys. 118, 2039 (2003)] with concepts generalizing the Born solution [Z. Phys. 1, 45 (1920)] for a solvated ion. First, we derive an exact representation according to which the sources of the RF potential and energy are inducible atomic anti-polarization densities and atomic shielding charge distributions. Modeling these atomic densities by Gaussians leads to an approximate representation. Here, the strengths of the Gaussian shielding charge distributions are directly given in terms of the static partial charges as defined, e.g., by standard MM force fields for the various atom types, whereas the strengths of the Gaussian anti-polarization densities are calculated by a self-consistency iteration. The atomic volumes are also described by Gaussians. To account for covalently overlapping atoms, their effective volumes are calculated by another self-consistency procedure, which guarantees that the dielectric function ε(r) is close to one everywhere inside the protein. The Gaussian widths σ(i) of the atoms i are parameters of the RF approximation. The remarkable accuracy of the method is demonstrated by comparison with Kirkwood's analytical solution for a spherical protein [J. Chem. Phys. 2, 351 (1934)] and with computationally expensive grid-based numerical solutions for simple model systems in dielectric continua including a di-peptide (Ac-Ala-NHMe) as modeled by a standard MM force field. The latter example shows how weakly the RF conformational free energy landscape depends on the parameters σ(i). A summarizing discussion highlights the achievements of the new theory and of its approximate solution particularly by comparison with so-called generalized Born methods. A follow-up paper describes how the method enables Hamiltonian, efficient, and accurate MM molecular dynamics simulations of proteins in dielectric solvent continua.

Electrostatics of proteins in dielectric solvent continua. II. Hamiltonian reaction field dynamics

In Paper I of this work [S. Bauer, G. Mathias, and P. Tavan, J. Chem. Phys. 140, 104102 (2014)] we have presented a reaction field (RF) method, which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of polarizable molecular mechanics (MM) force fields. Building upon these results, here we suggest a method for linearly scaling Hamiltonian RF/MM molecular dynamics (MD) simulations, which we call "Hamiltonian dielectric solvent" (HADES). First, we derive analytical expressions for the RF forces acting on the solute atoms. These forces properly account for all those conditions, which have to be self-consistently fulfilled by RF quantities introduced in Paper I. Next we provide details on the implementation, i.e., we show how our RF approach is combined with a fast multipole method and how the self-consistency iterations are accelerated by the use of the so-called direct inversion in the iterative subspace. Finally we demonstrate that the method and its implementation enable Hamiltonian, i.e., energy and momentum conserving HADES-MD, and compare in a sample application on Ac-Ala-NHMe the HADES-MD free energy landscape at 300 K with that obtained in Paper I by scanning of configurations and with one obtained from an explicit solvent simulation.

Understanding nucleic acid-ion interactions

Ions surround nucleic acids in what is referred to as an ion atmosphere. As a result, the folding and dynamics of RNA and DNA and their complexes with proteins and with each other cannot be understood without a reasonably sophisticated appreciation of these ions' electrostatic interactions. However, the underlying behavior of the ion atmosphere follows physical rules that are distinct from the rules of site binding that biochemists are most familiar and comfortable with. The main goal of this review is to familiarize nucleic acid experimentalists with the physical concepts that underlie nucleic acid-ion interactions. Throughout, we provide practical strategies for interpreting and analyzing nucleic acid experiments that avoid pitfalls from oversimplified or incorrect models. We briefly review the status of theories that predict or simulate nucleic acid-ion interactions and experiments that test these theories. Finally, we describe opportunities for going beyond phenomenological fits to a next-generation, truly predictive understanding of nucleic acid-ion interactions.

Molecular modeling study on the allosteric inhibition mechanism of HIV-1 integrase by LEDGF/p75 binding site inhibitors

HIV-1 integrase (IN) is essential for the integration of viral DNA into the host genome and an attractive therapeutic target for developing antiretroviral inhibitors. LEDGINs are a class of allosteric inhibitors targeting LEDGF/p75 binding site of HIV-1 IN. Yet, the detailed binding mode and allosteric inhibition mechanism of LEDGINs to HIV-1 IN is only partially understood, which hinders the structure-based design of more potent anti-HIV agents. A molecular modeling study combining molecular docking, molecular dynamics simulation, and binding free energy calculation were performed to investigate the interaction details of HIV-1 IN catalytic core domain (CCD) with two recently discovered LEDGINs BI-1001 and CX14442, as well as the LEDGF/p75 protein. Simulation results demonstrated the hydrophobic domain of BI-1001 and CX14442 engages one subunit of HIV-1 IN CCD dimer through hydrophobic interactions, and the hydrophilic group forms hydrogen bonds with HIV-1 IN CCD residues from other subunit. CX14442 has a larger tert-butyl group than the methyl of BI-1001, and forms better interactions with the highly hydrophobic binding pocket of HIV-1 IN CCD dimer interface, which can explain the stronger affinity of CX14442 than BI-1001. Analysis of the binding mode of LEDGF/p75 with HIV-1 IN CCD reveals that the LEDGF/p75 integrase binding domain residues Ile365, Asp366, Phe406 and Val408 have significant contributions to the binding of the LEDGF/p75 to HIV1-IN. Remarkably, we found that binding of BI-1001 and CX14442 to HIV-1 IN CCD induced the structural rearrangements of the 140 s loop and oration displacements of the side chains of the three conserved catalytic residues Asp64, Asp116, and Glu152 located at the active site. These results we obtained will be valuable not only for understanding the allosteric inhibition mechanism of LEDGINs but also for the rational design of allosteric inhibitors of HIV-1 IN targeting LEDGF/p75 binding site.

Computational study on the inhibitor binding mode and allosteric regulation mechanism in hepatitis C virus NS3/4A protein

HCV NS3/4A protein is an attractive therapeutic target responsible for harboring serine protease and RNA helicase activities during the viral replication. Small molecules binding at the interface between the protease and helicase domains can stabilize the closed conformation of the protein and thus block the catalytic function of HCV NS3/4A protein via an allosteric regulation mechanism. But the detailed mechanism remains elusive. Here, we aimed to provide some insight into the inhibitor binding mode and allosteric regulation mechanism of HCV NS3/4A protein by using computational methods. Four simulation systems were investigated. They include: apo state of HCV NS3/4A protein, HCV NS3/4A protein in complex with an allosteric inhibitor and the truncated form of the above two systems. The molecular dynamics simulation results indicate HCV NS3/4A protein in complex with the allosteric inhibitor 4VA adopts a closed conformation (inactive state), while the truncated apo protein adopts an open conformation (active state). Further residue interaction network analysis suggests the communication of the domain-domain interface play an important role in the transition from closed to open conformation of HCV NS3/4A protein. However, the inhibitor stabilizes the closed conformation through interaction with several key residues from both the protease and helicase domains, including His57, Asp79, Asp81, Asp168, Met485, Cys525 and Asp527, which blocks the information communication between the functional domains interface. Finally, a dynamic model about the allosteric regulation and conformational changes of HCV NS3/4A protein was proposed and could provide fundamental insights into the allosteric mechanism of HCV NS3/4A protein function regulation and design of new potent inhibitors.

Molecular modeling and residue interaction network studies on the mechanism of binding and resistance of the HCV NS5B polymerase mutants to VX-222 and ANA598

Hepatitis C virus (HCV) NS5B protein is an RNA-dependent RNA polymerase (RdRp) with essential functions in viral genome replication and represents a promising therapeutic target to develop direct-acting antivirals (DAAs). Multiple nonnucleoside inhibitors (NNIs) binding sites have been identified within the polymerase. VX-222 and ANA598 are two NNIs targeting thumb II site and palm I site of HCV NS5B polymerase, respectively. These two molecules have been shown to be very effective in phase II clinical trials. However, the emergence of resistant HCV replicon variants (L419M, M423T, I482L mutants to VX-222 and M414T, M414L, G554D mutants to ANA598) has significantly decreased their efficacy. To elucidate the molecular mechanism about how these mutations influenced the drug binding mode and decreased drug efficacy, we studied the binding modes of VX-222 and ANA598 to wild-type and mutant polymerase by molecular modeling approach. Molecular dynamics (MD) simulations results combined with binding free energy calculations indicated that the mutations significantly altered the binding free energy and the interaction for the drugs to polymerase. The further per-residue binding free energy decomposition analysis revealed that the mutations decreased the interactions with several key residues, such as L419, M423, L474, S476, I482, L497, for VX-222 and L384, N411, M414, Y415, Q446, S556, G557 for ANA598. These were the major origins for the resistance to these two drugs. In addition, by analyzing the residue interaction network (RIN) of the complexes between the drugs with wild-type and the mutant polymerase, we found that the mutation residues in the networks involved in the drug resistance possessed a relatively lower size of topology centralities. The shift of betweenness and closeness values of binding site residues in the mutant polymerase is relevant to the mechanism of drug resistance of VX-222 and ANA598. These results can provide an atomic-level understanding about the mechanisms of drug resistance conferred by the studied mutations and will be helpful to design more potent inhibitors which could effectively overcome drug resistance of antivirus agents.

Inhibition of antithrombin and bovine serum albumin native state aggregation by heparin

Protein native state aggregation, a major problem in pharmaceutical and biological processes, has been addressed pharmacologically by the addition of protein-binding excipients. Heparin (Hp), a highly sulfated polysaccharide, interacts with numerous proteins with moderate to high affinity, but reports about its effect on protein aggregation are contradictory. We studied the pH dependence of the aggregation of antithrombin (AT) and bovine serum albumin (BSA) in the presence and absence of heparin. High-precision turbidimetry showed strong aggregation for both AT and BSA in I = 10 mM NaCl, conditions at which electrostatically driven Hp binding and aggregation both occur, with more obvious aggregation of heparin-free AT appearing as larger aggregate size. Aggregation of AT was dramatically inhibited at Hp: protein 6:1 (mole ratio); however, the effect at 0.5:1 Hp:protein was greater for BSA. Frontal analysis capillary electrophoresis showed a much larger equilibrium association constant Kobs between Hp and AT, in accord with the onset of Hp binding at a higher pH; both effects are explained by the higher charge density of the positive domain for AT as revealed by modeling with DelPhi. The corresponding modeling images showed that these domains persist at high salt only for AT, consistent with the 160-fold drop in Kobs at 100 mM salt for BSA-Hp binding. The smaller inhibition effect for AT arises from the tendency of its uncomplexed monomer to form larger aggregates more rapidly, but the stronger binding of Hp to AT does not facilitate Hp-induced aggregate dissolution which occurs more readily for BSA. This can be attributed to the higher density of AT aggregates evidenced by higher fractal dimensions. Differences between inhibition and reversal by Hp arise because the former may depend on the stage at which Hp enters the aggregation process and the latter on aggregate size and morphology.

A Poisson-Boltzmann Equation Test Model for Protein in Spherical Solute Region and its Applications

The Poisson-Boltzmann equation (PBE) is one important implicit solvent continuum model for calculating electrostatics of protein in ionic solvent. Several numerical algorithms and program packages have been developed but verification and comparison between them remains an interesting topic. In this paper, a PBE test model is presented for a protein in a spherical solute region, along with its analytical solution. It is then used to verify a PBE finite element solver and applied to a numerical comparison study between a finite element solver and a finite difference solver. Such a study demonstrates the importance of retaining the interface conditions in the development of PBE solvers.

Redox specificity of 2-hydroxyacid-coupled NAD(+)/NADH dehydrogenases: a study exploiting "reactive" arginine as a reporter of protein electrostatics

With "reactive" arginine as a kinetic reporter, 2-hydroxyacid dehydrogenases are assessed in basis of their specialization as NAD(+)-reducing or NADH-oxidizing enzymes. Specifically, M4 and H4 lactate dehydrogenases (LDHs) and cytoplasmic and mitochondrial malate dehydrogenases (MDHs) are compared to assess if their coenzyme specificity may involve electrostatics of cationic or neutral nicotinamide structure as the basis. The enzymes from diverse eukaryote and prokaryote sources thus are assessed in "reactivity" of functionally-critical arginine as a function of salt concentration and pH. Electrostatic calculations were performed on "reactive" arginines and found good correspondence with experiment. The reductive and oxidative LDHs and MDHs are assessed in their count over ionizable residues and in placement details of the residues in their structures as proteins. The variants found to be high or low in ΔpKa of "reactive" arginine are found to be also strong or weak cations that preferentially oxidize NADH (neutral nicotinamide structure) or reduce NAD(+) (cationic nicotinamide structure). The ionized groups of protein structure may thus be important to redox specificity of the enzyme on basis of electrostatic preference for the oxidized (cationic nicotinamide) or reduced (neutral nicotinamide) coenzyme. Detailed comparisons of isozymes establish that the residues contributing in their redox specificity are scrambled in structure of the reductive enzyme.

Multiscale Multiphysics and Multidomain Models I: Basic Theory

This work extends our earlier two-domain formulation of a differential geometry based multiscale paradigm into a multidomain theory, which endows us the ability to simultaneously accommodate multiphysical descriptions of aqueous chemical, physical and biological systems, such as fuel cells, solar cells, nanofluidics, ion channels, viruses, RNA polymerases, molecular motors and large macromolecular complexes. The essential idea is to make use of the differential geometry theory of surfaces as a natural means to geometrically separate the macroscopic domain of solvent from the microscopic domain of solute, and dynamically couple continuum and discrete descriptions. Our main strategy is to construct energy functionals to put on an equal footing of multiphysics, including polar (i.e., electrostatic) solvation, nonpolar solvation, chemical potential, quantum mechanics, fluid mechanics, molecular mechanics, coarse grained dynamics and elastic dynamics. The variational principle is applied to the energy functionals to derive desirable governing equations, such as multidomain Laplace-Beltrami (LB) equations for macromolecular morphologies, multidomain Poisson-Boltzmann (PB) equation or Poisson equation for electrostatic potential, generalized Nernst-Planck (NP) equations for the dynamics of charged solvent species, generalized Navier-Stokes (NS) equation for fluid dynamics, generalized Newton's equations for molecular dynamics (MD) or coarse-grained dynamics and equation of motion for elastic dynamics. Unlike the classical PB equation, our PB equation is an integral-differential equation due to solvent-solute interactions. To illustrate the proposed formalism, we have explicitly constructed three models, a multidomain solvation model, a multidomain charge transport model and a multidomain chemo-electro-fluid-MD-elastic model. Each solute domain is equipped with distinct surface tension, pressure, dielectric function, and charge density distribution. In addition to long-range Coulombic interactions, various non-electrostatic solvent-solute interactions are considered in the present modeling. We demonstrate the consistency between the non-equilibrium charge transport model and the equilibrium solvation model by showing the systematical reduction of the former to the latter at equilibrium. This paper also offers a brief review of the field.

Structure of the uracil complex of Vaccinia virus uracil DNA glycosylase

Poxvirus uracil DNA glycosylases are the most diverse members of the family I uracil DNA glycosylases (UNGs). The crystal structure of the uracil complex of Vaccinia virus uracil DNA glycosylase (D4) was determined at 2.03 Å resolution. One uracil molecule was located in the active-site pocket in each of the 12 noncrystallographic symmetry-related D4 subunits. Although the UNGs of the poxviruses (including D4) feature significant differences in the characteristic motifs designated for uracil recognition and in the base-excision mechanism, the architecture of the active-site pocket in D4 is very similar to that in UNGs of other organisms. Overall, the interactions of the bound uracil with the active-site residues are also similar to the interactions previously observed in the structures of human and Escherichia coli UNG.

GPU linear and non-linear Poisson-Boltzmann solver module for DelPhi

SUMMARY

In this work, we present a CUDA-based GPU implementation of a Poisson-Boltzmann equation solver, in both the linear and non-linear versions, using double precision. A finite difference scheme is adopted and made suitable for the GPU architecture. The resulting code was interfaced with the electrostatics software for biomolecules DelPhi, which is widely used in the computational biology community. The algorithm has been implemented using CUDA and tested over a few representative cases of biological interest. Details of the implementation and performance test results are illustrated. A speedup of ~10 times was achieved both in the linear and non-linear cases.

AVAILABILITY AND IMPLEMENTATION

The module is open-source and available at http://www.electrostaticszone.eu/index.php/downloads.

Structure of the Toxoplasma gondii ROP18 kinase domain reveals a second ligand binding pocket required for acute virulence

At least a third of the human population is infected with the intracellular parasite Toxoplasma gondii, which contributes significantly to the disease burden in immunocompromised and neutropenic hosts and causes serious congenital complications when vertically transmitted to the fetus. Genetic analyses have identified the Toxoplasma ROP18 Ser/Thr protein kinase as a major factor mediating acute virulence in mice. ROP18 is secreted into the host cell during the invasion process, and its catalytic activity is required for the acute virulence phenotype. However, its precise molecular function and regulation are not fully understood. We have determined the crystal structure of the ROP18 kinase domain, which is inconsistent with a previously proposed autoinhibitory mechanism of regulation. Furthermore, a sucrose molecule bound to our structure identifies an additional ligand-binding pocket outside of the active site cleft. Mutational analysis confirms an important role for this pocket in virulence.

Exploring accurate Poisson-Boltzmann methods for biomolecular simulations

Accurate and efficient treatment of electrostatics is a crucial step in computational analyses of biomolecular structures and dynamics. In this study, we have explored a second-order finite-difference numerical method to solve the widely used Poisson-Boltzmann equation for electrostatic analyses of realistic bio-molecules. The so-called immersed interface method was first validated and found to be consistent with the classical weighted harmonic averaging method for a diversified set of test biomolecules. The numerical accuracy and convergence behaviors of the new method were next analyzed in its computation of numerical reaction field grid potentials, energies, and atomic solvation forces. Overall similar convergence behaviors were observed as those by the classical method. Interestingly, the new method was found to deliver more accurate and better-converged grid potentials than the classical method on or nearby the molecular surface, though the numerical advantage of the new method is reduced when grid potentials are extrapolated to the molecular surface. Our exploratory study indicates the need for further improving interpolation/extrapolation schemes in addition to the developments of higher-order numerical methods that have attracted most attention in the field.

Heparin decamer bridges a growth factor and an oligolysine by different charge-driven interactions

Full-length heparin is widely used in tissue engineering applications due its multiple protein-binding sites that allow it to retain growth factor affinity while associating with oligopeptide components of the tissue scaffold. However, the extent to which oligopeptide coupling interferes with cognate protein binding is difficult to predict. To investigate such simultaneous interactions, we examined a well-defined ternary system comprised of acidic fibroblast growth factor (FGF), tetralysine (K4), with a heparin decamer (dp10) acting as a noncovalent coupler. Electrospray ionization mass spectrometry was used to assess binding affinities and complex stoichiometries as a function of ionic strength for dp10·K4 and FGF·dp10. The ionic strength dependence of K4·dp10 formation is qualitatively consistent with binding driven by the release of condensed counterions previously suggested for native heparin with divalent oligopeptides (Mascotti, D. P.; Lohman, T. M. Biochemistry 1995, 34, 2908-2915). On the other hand, FGF binding displays more complex ionic strength dependence, with higher salt resistance. Remarkably, dp10 that can bind two FGF molecules can only bind one tetralysine. The limited binding of K4 to dp10 suggests that the tetralysine might not block growth factor binding, and the 1:1:1 ternary complex is indeed observed. The analysis of mass distribution of the bound dp10 chains in FGF·dp10, FGF2·dp10, and FGF·dp10·K4 complexes indicated that higher degrees of dp10 sulfation promote the formation of FGF2·dp10 and FGF·dp10·K4. Thus, the selectivity of appropriately chosen short heparin chains could be used to modulate growth factor sequestration and release in a way not feasible with heterogeneous native heparin. In support of this, human hepatocellular carcinoma cells (HEP3Bs) treated with FGF·dp10·K4 were found to exhibit biological activity similar to cells treated with FGF.

Binding selectivity studies of PKBα using molecular dynamics simulation and free energy calculations

Designing selective protein kinase B (PKB/Akt) inhibitor is an area of intense research to develop potential anticancer drugs. In the present study, the molecular basis governing PKB-selective inhibition has been investigated using molecular dynamics simulation. The binding free energies calculated by MM/PBSA gave a good correlation with the experimental biological activity and a good explanation of the activity difference of the studied inhibitors. The decomposition of free energies by MM/GBSA indicates that the ethyl group on pyrrolo[2,3-d]pyrimidine ring of inhibitor Lig1 (N-{[(3S)-3-amino-1-(5-ethyl-7H-pyrrolo[2,3-d]pyrimidin-4-yl)pyrrolidin-3-yl]-methyl}-2,4-difluoro-benzamide) is an important contributor to its PKBα selectivity due to its hydrophobic interaction with the side chain of Thr291 in PKBα. The substituted groups on the pyrrolidine ring of Lig1 also show a strong tendency to mediate protein-ligand interactions through the hydrogen bonds formed between the amino or amide groups of Lig1 and the carboxyl O atoms of Glu234, Glu278, and Asp292 of PKBα. It was reported that there are only three key amino acid differences between PKBα (Thr211, Ala230, Met281) and PKA (Val104, Val123, Leu173) within the clefts of ATP-binding sites. These differences propel a drastic conformational change in PKA, weakening its binding interactions with inhibitor. The impact was also confirmed by MD simulated interaction modes of inhibitor binding to PKBα mutants with the in silico mutations of the three key amino acids, respectively. We expect that the results obtained here could be useful for future rational design of specific ATP-competitive inhibitors of PKBα.

Engineering of weak helper interactions for high-efficiency FRET probes

Fluorescence resonance energy transfer (FRET)-based detection of protein interactions is limited by the very narrow range of FRET-permitting distances. We show two different strategies for the rational design of weak helper interactions that co-recruit donor and acceptor fluorophores for a more robust detection of bimolecular FRET: (i) in silico design of electrostatically driven encounter complexes and (ii) fusion of tunable domain-peptide interaction modules based on WW or SH3 domains. We tested each strategy for optimization of FRET between (m)Citrine and mCherry, which do not natively interact. Both approaches yielded comparable and large increases in FRET efficiencies with little or no background. Helper-interaction modules can be fused to any pair of fluorescent proteins and could, we found, enhance FRET between mTFP1 and mCherry as well as between mTurquoise2 and mCitrine. We applied enhanced helper-interaction FRET (hiFRET) probes to study the binding between full-length H-Ras and Raf1 as well as the drug-induced interaction between Raf1 and B-Raf.

Quantum mechanical modeling: a tool for the understanding of enzyme reactions

Most enzyme reactions involve formation and cleavage of covalent bonds, while electrostatic effects, as well as dynamics of the active site and surrounding protein regions, may also be crucial. Accordingly, special computational methods are needed to provide an adequate description, which combine quantum mechanics for the reactive region with molecular mechanics and molecular dynamics describing the environment and dynamic effects, respectively. In this review we intend to give an overview to non-specialists on various enzyme models as well as established computational methods and describe applications to some specific cases. For the treatment of various enzyme mechanisms, special approaches are often needed to obtain results, which adequately refer to experimental data. As a result of the spectacular progress in the last two decades, most enzyme reactions can be quite precisely treated by various computational methods.

Isolating influential regions of electrostatic focusing in protein and DNA structure

Electrostatic focusing is a general phenomenon that occurs in cavities and grooves on the molecular surface of biomolecules. Narrow surface features can partially shield charged atoms from the high-dielectric solvent, enhancing electrostatic potentials inside the cavity and projecting electric field lines outward into the solvent. This effect has been observed in many instances and is widely considered in the human examination of molecular structure, but it is rarely integrated into the digital representations used in protein structure comparison software. To create a computational representation of electrostatic focusing, that is compatible with structure comparison algorithms, this paper presents an approach that generates three-dimensional solids that approximate regions where focusing occurs. We verify the accuracy of this representation against instances of focusing in proteins and DNA. Noting that this representation also identifies thin focusing regions on the molecular surface that are unlikely to affect binding, we describe a second algorithm that conservatively isolates larger focusing regions. The resulting 3D solids can be compared with Boolean set operations, permitting a new range of analyses on the regions where electrostatic focusing occurs. They also represent a novel integration of molecular shape and electrostatic focusing into the same structure comparison framework.

Computational Alanine Scanning Mutagenesis-An Improved Methodological Approach for Protein-DNA Complexes

Proteins and protein-based complexes are the basis of many key systems in nature and have been the subject of intense research in the last decades, in an attempt to acquire comprehensive knowledge of reactions that take place in nature. Computational Alanine Scanning Mutagenesis approaches have been extensively used in the study of protein interfaces and in the determination of the most important residues for complex formation, the Hot-spots. However, as it is usually applied to the study of protein-protein interfaces, we tried to modify and apply it to the study of protein-DNA interfaces, which are also crucial in nature but have not been the subject of as much research. In this work, we carry out MD simulations of seven protein-DNA complexes and tested the influence of the variation of different parameters on the determination of the binding free energy terms (ΔΔGbinding) of 78 mutations: solvent representation, internal dielectric constant, Linear and Nonlinear Poisson-Boltzmann equation, Generalized Born model, simulation time, number of structures analyzed, number of MD trajectories, force field used, and energetic terms involved. Overall, this new approach gave an average error of 1.55 kcal/mol, and P, R, F1, accuracy, and specificity values of 0.78, 0.50, 0.61, 0.77, and 0.92, respectively. This improved computational alanine scanning mutagenesis approach may serve as a tool to explore the behavior of this important class of complexes.

Distance-independent charge recombination kinetics in cytochrome c-cytochrome c peroxidase complexes: compensating changes in the electronic coupling and reorganization energies

Charge recombination rate constants vary no more than 3-fold for interprotein ET in the Zn-substituted wild type (WT) cytochrome c peroxidase (CcP):cytochrome c (Cc) complex and in complexes with four mutants of the Cc protein (i.e., F82S, F82W, F82Y, and F82I), despite large differences in the ET distance. Theoretical analysis indicates that charge recombination for all complexes involves a combination of tunneling and hopping via Trp191. For three of the five structures (WT and F82S(W)), the protein favors hopping more than that in the other two structures that have longer heme → ZnP distances (F82Y(I)). Experimentally observed biexponential ET kinetics is explained by the complex locking in alternative coupling pathways, where the acceptor hole state is either primarily localized on ZnP (slow phase) or on Trp191 (fast phase). The large conformational differences between the CcP:Cc interface for the F82Y(I) mutants compared to that the WT and F82S(W) complexes are predicted to change the reorganization energies for the CcP:Cc ET reactions because of changes in solvent exposure and interprotein ET distances. Since the recombination reaction is likely to occur in the inverted Marcus regime, an increased reorganization energy compensates the decreased role for hopping recombination (and the longer transfer distance) in the F82Y(I) mutants. Taken together, coupling pathway and reorganization energy effects for the five protein complexes explain the observed insensitivity of recombination kinetics to donor-acceptor distance and docking pose and also reveals how hopping through aromatic residues can accelerate long-range ET.

Molecular dynamics studies on the conformational transitions of adenylate kinase: a computational evidence for the conformational selection mechanism

Escherichia coli adenylate kinase (ADK) is a monomeric phosphotransferase enzyme that catalyzes reversible transfer of phosphoryl group from ATP to AMP with a large-scale domain motion. The detailed mechanism for this conformational transition remains unknown. In the current study, we performed long time-scale molecular dynamics simulations on both open and closed states of ADK. Based on the structural analyses of the simulation trajectories, we detected over 20 times conformational transitions between the open and closed states of ADK and identified two novel conformations as intermediate states in the catalytic processes. With these findings, we proposed a possible mechanism for the large-scale domain motion of Escherichia coli ADK and its catalytic process: (1) the substrate free ADK adopted an open conformation; (2) ATP bound with LID domain closure; (3) AMP bound with NMP domain closure; (4) phosphoryl transfer occurred with ATP, and AMP converted into two ADPs, and no conformational transition was detected in the enzyme; (5) LID domain opened with one ADP released; (6) another ADP released with NMP domain open. As both open and closed states sampled a wide range of conformation transitions, our simulation strongly supported the conformational selection mechanism for Escherichia coli ADK.

Structural characterization of a unique marine animal family 7 cellobiohydrolase suggests a mechanism of cellulase salt tolerance

Nature uses a diversity of glycoside hydrolase (GH) enzymes to convert polysaccharides to sugars. As lignocellulosic biomass deconstruction for biofuel production remains costly, natural GH diversity offers a starting point for developing industrial enzymes, and fungal GH family 7 (GH7) cellobiohydrolases, in particular, provide significant hydrolytic potential in industrial mixtures. Recently, GH7 enzymes have been found in other kingdoms of life besides fungi, including in animals and protists. Here, we describe the in vivo spatial expression distribution, properties, and structure of a unique endogenous GH7 cellulase from an animal, the marine wood borer Limnoria quadripunctata (LqCel7B). RT-quantitative PCR and Western blot studies show that LqCel7B is expressed in the hepatopancreas and secreted into the gut for wood degradation. We produced recombinant LqCel7B, with which we demonstrate that LqCel7B is a cellobiohydrolase and obtained four high-resolution crystal structures. Based on a crystallographic and computational comparison of LqCel7B to the well-characterized Hypocrea jecorina GH7 cellobiohydrolase, LqCel7B exhibits an extended substrate-binding motif at the tunnel entrance, which may aid in substrate acquisition and processivity. Interestingly, LqCel7B exhibits striking surface charges relative to fungal GH7 enzymes, which likely results from evolution in marine environments. We demonstrate that LqCel7B stability and activity remain unchanged, or increase at high salt concentration, and that the L. quadripunctata GH mixture generally contains cellulolytic enzymes with highly acidic surface charge compared with enzymes derived from terrestrial microbes. Overall, this study suggests that marine cellulases offer significant potential for utilization in high-solids industrial biomass conversion processes.

Binding free-energy calculation of an ion-peptide complex by constrained dynamics

Binding free energy is the most important physical parameter that describes the binding affinity of a receptor-ligand complex. Conventionally, it was obtained based on the thermodynamic cycle or alchemical reaction. These strategies have been widely used, but they would be problematic if the receptors and/or ligands have large conformational changes during the binding processes. In this paper, we present a way to calculate the binding free energy: constrained dynamics along a fragmental and high-dimensional transition path. This method directly considers unbound states in the simulation. The application to the calmodulin loop-calcium complexes shows that it is practical and the calculated relative binding affinities are in good agreement with experimental results.