Benzhuo Lu

Computational analysis and prediction of the binding motif and protein interacting partners of the Abl SH3 domain

Protein-protein interactions, particularly weak and transient ones, are often mediated by peptide recognition domains, such as Src Homology 2 and 3 (SH2 and SH3) domains, which bind to specific sequence and structural motifs. It is important but challenging to determine the binding specificity of these domains accurately and to predict their physiological interacting partners. In this study, the interactions between 35 peptide ligands (15 binders and 20 non-binders) and the Abl SH3 domain were analyzed using molecular dynamics simulation and the Molecular Mechanics/Poisson-Boltzmann Solvent Area method. The calculated binding free energies correlated well with the rank order of the binding peptides and clearly distinguished binders from non-binders. Free energy component analysis revealed that the van der Waals interactions dictate the binding strength of peptides, whereas the binding specificity is determined by the electrostatic interaction and the polar contribution of desolvation. The binding motif of the Abl SH3 domain was then determined by a virtual mutagenesis method, which mutates the residue at each position of the template peptide relative to all other 19 amino acids and calculates the binding free energy difference between the template and the mutated peptides using the Molecular Mechanics/Poisson-Boltzmann Solvent Area method. A single position mutation free energy profile was thus established and used as a scoring matrix to search peptides recognized by the Abl SH3 domain in the human genome. Our approach successfully picked ten out of 13 experimentally determined binding partners of the Abl SH3 domain among the top 600 candidates from the 218,540 decapeptides with the PXXP motif in the SWISS-PROT database. We expect that this physical-principle based method can be applied to other protein domains as well.

Order N algorithm for computation of electrostatic interactions in biomolecular systems

Poisson–Boltzmann electrostatics is a well established model in biophysics; however, its application to large-scale biomolecular processes such as protein–protein encounter is still limited by the efficiency and memory constraints of existing numerical techniques. In this article, we present an efficient and accurate scheme that incorporates recently developed numerical techniques to enhance our computational ability. In particular, a boundary integral equation approach is applied to discretize the linearized Poisson–Boltzmann equation; the resulting integral formulas are well conditioned and are extended to systems with arbitrary numbers of biomolecules. The solution process is accelerated by Krylov subspace methods and a new version of the fast multipole method. In addition to the electrostatic energy, fast calculations of the forces and torques are made possible by using an interpolation procedure. Numerical experiments show that the implemented algorithm is asymptotically optimal O(N) in both CPU time and required memory, and application to the acetylcholinesterase–fasciculin complex is illustrated.

Computation of electrostatic forces between solvated molecules determined by the Poisson–Boltzmann equation using a boundary element method

A rigorous approach is proposed to calculate the electrostatic forces among an arbitrary number of solvated molecules in ionic solution determined by the linearized Poisson-Boltzmann equation. The variational principle is used and implemented in the frame of a boundary element method (BEM). This approach does not require the calculation of the Maxwell stress tensor on the molecular surface, therefore it totally avoids the hypersingularity problem in the direct BEM whenever one needs to calculate the gradient of the surface potential or the stress tensor. This method provides an accurate and efficient way to calculate the full intermolecular electrostatic interaction energy and force, which could potentially be used in Brownian dynamics simulation of biomolecular association. The method has been tested on some simple cases to demonstrate its reliability and efficiency, and parts of the results are compared with analytical results and with those obtained by some known methods such as adaptive Poisson-Boltzmann solver.

Electrodiffusion: A continuum modeling framework for biomolecular systems with realistic spatiotemporal resolution

A computational framework is presented for the continuum modeling of cellular biomolecular diffusion influenced by electrostatic driving forces. This framework is developed from a combination of state-of-the-art numerical methods, geometric meshing, and computer visualization tools. In particular, a hybrid of (adaptive) finite element and boundary element methods is adopted to solve the Smoluchowski equation (SE), the Poisson equation (PE), and the Poisson-Nernst-Planck equation (PNPE) in order to describe electrodiffusion processes. The finite element method is used because of its flexibility in modeling irregular geometries and complex boundary conditions. The boundary element method is used due to the convenience of treating the singularities in the source charge distribution and its accurate solution to electrostatic problems on molecular boundaries. Nonsteady-state diffusion can be studied using this framework, with the electric field computed using the densities of charged small molecules and mobile ions in the solvent. A solution for mesh generation for biomolecular systems is supplied, which is an essential component for the finite element and boundary element computations. The uncoupled Smoluchowski equation and Poisson-Boltzmann equation are considered as special cases of the PNPE in the numerical algorithm, and therefore can be solved in this framework as well. Two types of computations are reported in the results: stationary PNPE and time-dependent SE or Nernst-Planck equations solutions. A biological application of the first type is the ionic density distribution around a fragment of DNA determined by the equilibrium PNPE. The stationary PNPE with nonzero flux is also studied for a simple model system, and leads to an observation that the interference on electrostatic field of the substrate charges strongly affects the reaction rate coefficient. The second is a time-dependent diffusion process: the consumption of the neurotransmitter acetylcholine by acetylcholinesterase, determined by the SE and a single uncoupled solution of the Poisson-Boltzmann equation. The electrostatic effects, counterion compensation, spatiotemporal distribution, and diffusion-controlled reaction kinetics are analyzed and different methods are compared.

Poisson-Nernst-Planck Equations for Simulating Biomolecular Diffusion-Reaction Processes II: Size Effects on Ionic Distributions and Diffusion-Reaction Rates

The effects of finite particle size on electrostatics, density profiles, and diffusion have been a long existing topic in the study of ionic solution. The previous size-modified Poisson-Boltzmann and Poisson-Nernst-Planck models are revisited in this article. In contrast to many previous works that can only treat particle species with a single uniform size or two sizes, we generalize the Borukhov model to obtain a size-modified Poisson-Nernst-Planck (SMPNP) model that is able to treat nonuniform particle sizes. The numerical tractability of the model is demonstrated as well. The main contributions of this study are as follows. 1), We show that an (arbitrarily) size-modified PB model is indeed implied by the SMPNP equations under certain boundary/interface conditions, and can be reproduced through numerical solutions of the SMPNP. 2), The size effects in the SMPNP effectively reduce the densities of highly concentrated counterions around the biomolecule. 3), The SMPNP is applied to the diffusion-reaction process for the first time, to our knowledge. In the case of low substrate density near the enzyme reactive site, it is observed that the rate coefficients predicted by SMPNP model are considerably larger than those by the PNP model, suggesting both ions and substrates are subject to finite size effects. 4), An accurate finite element method and a convergent Gummel iteration are developed for the numerical solution of the completely coupled nonlinear system of SMPNP equations.

Free Energy for the Permeation of Na + and Cl - Ions and Their Ion-Pair through a Zwitterionic Dimyristoyl Phosphatidylcholine Lipid Bilayer by Umbrella Integration with Harmonic Fourier Beads

Understanding the mechanism of ion permeation across lipid bilayers is key to controlling osmotic pressure and developing new ways of delivering charged, drug-like molecules inside cells. Recent reports suggest ion-pairing as the mechanism to lower the free energy barrier for the ion permeation in disagreement with predictions from the simple electrostatic models. In this paper we quantify the effect of ion-pairing or charge quenching on the permeation of Na+ and Cl− ions across DMPC lipid bilayer by computing the corresponding potentials of mean force (PMFs) using fully atomistic molecular dynamics simulations. We find that the free energy barrier to permeation reduces in the order Na+−Cl− ion-pair (27.6 kcal/mol) > Cl− (23.6 kcal/mol) > Na+ (21.9 kcal/mol). Furthermore, with the help of these PMFs we derive the change in the binding free energy between the Na+ and Cl− with respect to that in water as a function of the bilayer permeation depth. Despite the fact that the bilayer boosts the Na+−Cl− ion binding free energy by as high as 17.9 kcal/mol near its center, ion-pairing between such hydrophilic ions as Na+ and Cl− does not assist their permeation. However, based on a simple thermodynamic cycle, we suggest that ion-pairing between ions of opposite charge and solvent philicity could enhance ion permeation. Comparison of the computed permeation barriers for Na+ and Cl− ions with available experimental data supports this notion. This work establishes general computational methodology to address ion-pairing in fluid anisotropic media and details the ion permeation mechanism on atomic level.

TMSmesh: A Robust Method for Molecular Surface Mesh Generation Using a Trace Technique.

Qualified, stable, and efficient molecular surface meshing appears to be necessitated by recent developments for realistic mathematical modeling and numerical simulation of biomolecules, especially in implicit solvent modeling (e.g., see a review in B. Z. Lu et al. Commun. Comput. Phys. 2008, 3, 973-1009). In this paper, we present a new method: tracing molecular surface for meshing (TMSmesh) the Gaussian surface of biomolecules. The method computes the surface points by solving a nonlinear equation directly, polygonizes by connecting surface points through a trace technique, and finally outputs a triangulated mesh. TMSmesh has a linear complexity with respect to the number of atoms and is shown to be capable of handling molecules consisting of more than one million atoms, which is usually difficult for the existing methods for surface generation used in molecular visualization and geometry analysis. Moreover, the meshes generated by TMSmesh are successfully tested in boundary element solutions of the Poisson-Boltzmann equation, which directly gives rise to a route to simulate electrostatic solvation of large-scale molecular systems. The binary version of TMSmesh and a set of representative PQR benchmark molecules are downloadable at our Web page http://lsec.cc.ac.cn/∼lubz/Meshing.html .

An Adaptive Fast Multipole Boundary Element Method for Poisson-Boltzmann Electrostatics

The numerical solution of the Poisson−Boltzmann (PB) equation is a useful but a computationally demanding tool for studying electrostatic solvation effects in chemical and biomolecular systems. Recently, we have described a boundary integral equation-based PB solver accelerated by a new version of the fast multipole method (FMM). The overall algorithm shows an order N complexity in both the computational cost and memory usage. Here, we present an updated version of the solver by using an adaptive FMM for accelerating the convolution type matrix-vector multiplications. The adaptive algorithm, when compared to our previous nonadaptive one, not only significantly improves the performance of the overall memory usage but also remarkably speeds the calculation because of an improved load balancing between the local- and far-field calculations. We have also implemented a node-patch discretization scheme that leads to a reduction of unknowns by a factor of 2 relative to the constant element method without sacrificing accuracy. As a result of these improvements, the new solver makes the PB calculation truly feasible for large-scale biomolecular systems such as a 30S ribosome molecule even on a typical 2008 desktop computer.

AFMPB: An adaptive fast multipole Poisson Boltzmann solver for calculating electrostatics in biomolecular systems

A Fortran program package is introduced for rapid evaluation of the electrostatic potentials and forces in biomolecular systems modeled by the linearized Poisson-Boltzmann equation. The numerical solver utilizes a well-conditioned boundary integral equation (BIE) formulation, a node-patch discretization scheme, a Krylov subspace iterative solver package with reverse communication protocols, and an adaptive new version of fast multipole method in which the exponential expansions are used to diagonalize the multipole to local translations. The program and its full description, as well as several closely related libraries and utility tools are available at http://lsec.cc.ac.cn/lubz/afmpb.html and a mirror site at http://mccammon.ucsd.edu/. This paper is a brief summary of the program: the algorithms, the implementation and the usage.

Protein molecular dynamics with electrostatic force entirely determined by a single Poisson-Boltzmann calculation.

The electrostatic force including the intramolecular Coulombic interactions and the electrostatic contribution of solvation effect were entirely calculated by using the finite difference Poisson‐Boltzmann method (FDPB), which was incorporated into the GROMOS96 force field to complete a new finite difference stochastic dynamics procedure (FDSD). Simulations were performed on an insulin dimer. Different relative dielectric constants were successively assigned to the protein interior; a value of 17 was selected as optimal for our system. The simulation data were analyzed and compared with those obtained from 500‐ps molecular dynamics (MD) simulation with explicit water and a 500‐ps conventional stochastic dynamics (SD) simulation without the mean solvent force. The results indicate that the FDSD method with GROMOS96 force field is suitable to study the dynamics and structure of proteins in solution if used with the optimal protein dielectric constant. Proteins 2002;48:497–504.