Understanding Protein Electrostatics Using Boundary-Integral Equations

1. Understanding Protein Electrostatics Using Boundary-Integral Equations Jaydeep P. Bardhan Dept. of Physiology and Molecular Biophysics Rush University Medical Center, Chicago IL Joint work with M. Knepley (Computation Institute, U. Chicago) P. Brune (Math and Computer Science Division, Argonne) A. Hildebrandt (J. Gutenberg U., Mainz, Germany)

2. Applied Math Biophysics My research Computer science (HPC) Outline: • Preliminaries: • Biomolecule electrostatics • Continuum theory and boundary-integral methods • Numerical simulation • Fast Poisson approximation • Nonlocal continuum model Emphasizing the interdisciplinary nature of computational biophysics

3. Fox L. Freberg Vander Kass ‘05 Fact: Water Makes Life Possible

4. d = 0 d = 1 A Crucial Consequence of Solvation • Molecular binding involves sacrificing solute--solvent interactions for solute--solute interactions:

5. Modeling ions in solution is critical! But today’s focus is on the simpler math of “pure” water. Basic Continuum Electrostatic Theory • 100-1000 times faster than MD • Protein model: • Shape: “union of spheres” (atoms) • Point charges at atom centers • Not very polarizable:  = 2-4 • Water model: no fixed charges • Single water: sphere of radius 1.4 Angstrom • Highly polarizable:  = 80 • In total: mixed-dielectric Poisson

6. PDE Complications Boundary conditions are at infinity Point charges must be spread onto the grid The dielectric interface is approximated Solving the PDE Directly is Possible, But… The idea: Just throw down a finite-difference grid or a finite-element mesh and go to town!

7. D Dirichlet: given Ex: Neumann: given (For exterior domain!) (3 dimensions) Potential anywhere in D Surface integrals ONLY Green’s Representation Formula • Well known: boundary values of a harmonic function determine it uniquely • The challenge is determining given BV. • Separation of variables • Numerical: finite elements, finite differences.. • In 3 dimensions, solve for 3-dimensional unknown • Alternatively: if you knew BOTH conditions, S Thus: finding the other boundary condition gives you the answer directly

8. Given data D S Deriving a Boundary Integral Equation: Exterior Neumann Problem • Given , need to find • Let r approach surface r is in the domain D r’ is on the boundary S

9. r r z D x + z S + y + x + + a + + + + y + + Continuous as + + + + + + + + + + + Also continuous as - - - - - - - - - - Addressing the Singularities • Single-layer potential • Dipole-layer Limit depends on WHICH SIDE of the surface your point is approaching!! a

10. D S Deriving a Boundary Integral Equation: Exterior Neumann Problem • Given , need to find

11. Medium problem: Exterior problems? Problems with mostly empty, uninteresting space To infinity Why Bother With Integral Equations? Easy problem:

12. PDE BIE Practical Advantages of PDE Approaches • Much more general (nonlinearity, etc) • Easier to parallelize (that’s different from “easy”) • Often easier to explore model space (see point 1) • PDE solvers give sparse systems; BIE, dense systems! • Less accessible (few codes) • Hard to convince it to run • Does one thing really well • Accessible (many codes) • Reliable, durable • Versatile, does OK job

13. 1 on panel i 0 elsewhere Enforce Enforce Galerkin: Similarity Between FEM and BEM • Both weighted residual methods: BEM FEM (Galerkin method)

14. Galerkin FEM: Smooth integrand: Easily computed with quadrature! Differences Between BEM and FEM • Extra freedom in choosing test functions • Matrix elements are harder to compute Collocation: test = delta functions Centroids of elements Galerkin BEM: Double integral of a singular function!!

15. Fast Solvers for Integral Equations Solve Ax=bapproximately using Krylov-subspace iterative methods such as GMRES: Compute dense matrix-vector product using O(N) method (fast multipole; tree code; precorrected FFT; FFTSVD) Improve iterative convergence with preconditioning For many problems, use diagonal entries! Iteration converges faster if matrix eigenvalues are “well clustered” P“looks like”A-1 Storing matrix: O(N2) time and memory Each multiply: O(N2) time Storing compressed matrix: O(N) time and memory Each multiply: O(N) time

16. + - + + + + - - + + - - + - - + Conservation law Constitutive relation A Boundary Integral Method For Biomolecule Electrostatics Boundary conditions handled exactly Point charges are treated exactly Meshing emphasis can be placed directly on the interface

17. BIBEE: A New, Rigorous Model of Continuum Electrostatics for Proteins “Boundary Integral Based Electrostatics Estimation” • Idea: Use preconditioner to approximate inverse • No need to compute sparsified operator (saves time and memory) • No need for Krylov solve • Test of elementary charges in a 20-Angstrom sphere: Single +1 charge +1, -1 charges 3 A apart

18. BIBEE: Introducing Different Variants The preconditioning approximation takes into account the singular character of the electric-field kernel: The Coulomb-field approximation ignores the operator entirely: CFA seems better here… …and worse here.

19.  i -1/10 -1/6 -1/2 A hundred years of analysis! BIBEE Approximates the Eigenvalues of the Boundary Integral Operator • The integral operator has to be split into two terms • BIBEE approximates E’s eigenvalues • P uses 0 (limit for sphere, prolate spheroid) • CFA uses -1/2 (known extremal) Sphere: analytical • Eigenvalues are real •  in [-1/2,+1/2) • -1/2 is always an EV • Left, right eigenvectors of -1/2 are constants

20. R1 R2 R3 BIBEE/CFA is the extension of CFA to multiple charges! No ad hoc parameters, no heuristic interpolation BIBEE Clarifies an Empirical, Heuristic Model + + BIBEE approx. charge includes all contributions “Effective Born radius” - the radius of a sphere with the same solvation energy Coulomb-field approximation: corresponds exactly to ignoring the integral operator. Still equation: the basis of totally nonphysical Generalized Born (GB) models Same approach taken by Borgis et al. in variational CFA

21. Feig et al. test set, > 600 proteins BIBEE/CFA Energy Is a Provable Upper Bound • BIBEE/P is an effective lower bound, provable in some cases but not all • Another variant (BIBEE/LB) is a provable LB but too loose to be useful Bardhan, Knepley, Anitescu (2009)

22.  i -1/10 -1/6 -1/2 BIBEE: Improve by Analyzing the Sphere • Get first mode (monopole) analytically correct, other modes are bounded from below: tighter lower bound! • Impact on sphere is better than impact on proteins (Feig et al. test set) Bardhan+Knepley, J. Chem. Phys. (in press)

23.  i -1/10 -1/6 Dominant energies come from dominant modes: try to capture dipole/quadrupole modes approximately! -1/2 BIBEE: Accurate One-parameter Model • This effective parameter is expected to be rigorously determined by approximating protein as ellipsoid (Onufriev+Sigalov, ‘06) Bardhan+Knepley, J. Chem. Phys. (in press)

24. Tripeptide Protein-Drug BIBEE: A New, Rigorous Model # boundary elements Total BEM time Matrix compression time SGB/CFA (heuristic) time BIBEE time • BIBEE is 3-5 times faster than full solve (including large setup time for both) • Unoptimized implementation (will save big on setup time) • Modern FMM implementation (Yokota, Knepley, Barba, et al.) gives 10-20X speedup

25. Reaction-Potential Operator Eigenvectors Have Physical Meaning • Eigenvectors from distinct eigenvalues are orthogonal • Thus: the eigenvectors correspond to charge distributions that do not interact via solvent polarization (weird, huh?) • If an approximate method generates a solvation matrix , its eigenvectors should “line up” well with the actual eigenvectors, i.e. i = j

26. “Getting the Modes Right” Is Important Modes from small eigenvalues still contribute significantly to the total energy Here, 25% of the total energy comes from modes with eigenvalues smaller than 1% of the maximum eigenvalue Eigenvalue Magnitude 104 -10 102 Projection of charge distribution onto eigenvector Cumulative Electrostatic Free Energy (kcal/mol) 100 -20 10-2 20 40 60 80 -30 Eigenvalue Index 20 40 60 80 Eigenvalue Index

27. SGB/CFA GBMV BIBEE/CFA BIBEE Is An Accurate, Parameter-Free Model Snapshots from MD • Peptide example Met-enkephalin BIBEE’s stronger “diagonal” appearance indicates superior reproduction of the eigenvectors of the operator. All models look essentially the same here.

28. Biophysics Computer science (HPC) BIBEE: A New, Rigorous Model of Continuum Electrostatics for Proteins Design systematic approximation Have proved that the model: • Gives upper and lower bounds • Preserves important physics Relate empirical models to strong math Leverages existing algorithms (e.g. fast multipole methods, parallel codes) Next: Apply to other physics problems Applied Math BIBEE

29. Lone pair electrons Hydrogen bonds Oxygen Hydrogens Nonlocal Continuum Electrostatics: Adding molecular realism “the right way” KNOWN weaknesses of Poisson model: • Linear response assumption Caveat: Nonlinearity IS important for more highly charged species! • Violates continuum-length-scale assumption • Water molecules have finite size • Water molecules form semi-structured networks First look for ways to extend existing models--don’t just give up and reinvent everything! Test with all-atom molecular dynamics Relatively small deviation! y=x denotes exactly linear response Nina, Beglov, Roux ‘97

30. Significant structuring of charge density! Data points: radii from molecular simulation (Aqvist 1990) and energies from experimental data Nonlocal Continuum Electrostatics: Demonstrating the Failure Mode Run all-atom molecular dynamics: ion surrounded by water Consequence: ion energies are wrong Ion radius in nanometers

31. Problems whose length scales are NOT well separated from those of the constituent molecules! Expect nonlocal theory to play major roles in nanoscale science and engineering modeling… de Abajo ‘08 Schatz et al. ‘01 Park ‘06 Duan et al. ‘07 Scott et al. ‘04 Gao et al., ‘09 Nonlocal Continuum Modeling:A Classical Multiscale Theory • Studied since the 1970s in numerous domains

32. Local response Wave number (inverse distance) Nonlocal Continuum Electrostatics:Nonlocal Dielectric Response • Polarization charge as a function of distance from the ion: not simple • Short-range: electronic response • Long-range: bulk behavior • Local: bulk everywhere • Nonlocal: simple function that captures asymptotes Supported by experiments and detailed simulations Smoothly interpolates between known limits

33. Green’s function for Nonlocal Continuum Electrostatics:Lorentzian Model and Promising Tests • Nonlocal response: • Now • Integrodifferential Poisson equation Single parameter fit for  gives much better agreement with experiment!!

34. Molecular surface “Licorice” “Cartoon” Nonlocal Continuum Electrostatics:Reformulation for Fast Simulations • Integrodifferential equations in complex geometries? • Result: No progress on nonlocal model for DECADES Spherical ions, charges near planar half-spaces… nothing else. • Breakthrough in 2004 (Hildebrandt et al.): • Define an auxiliary field: the displacement potential • Approximate the nonlocal boundary condition • Double reciprocity leads to a boundary-integral method

35. Yukawa/linearized Poisson-Boltzmann equation Displacement potential acts as a volume source Nonlocal Continuum Electrostatics:1. Introduce an Auxiliary Potential • Use Helmholtz decomposition: • Electrostatic potential now satisfies a Yukawa equation:

36. Nonlocal boundary condition: Choose to drop Nonlocal Continuum Electrostatics:2. Approximate Nonlocal B.C. • Original boundary conditions: • Exact normal deriv. of solvent potentials satisfy • The actual PDEs complete the local formulation:

37. Nonlocal Continuum Electrostatics:3. Green’s Theorem + Double Reciprocity • Electric potential Green’s theorem gives a volume integral • The displacement potential is harmonic: • Defining single- and double-layer operators 0

38. Nonlocal Continuum Electrostatics: Purely BIE Formulation • Three surface variables, two types of Green’s functions, and a mixed first-second kind problem • Fasel et al. have recently derived a purely second-kind method Hildebrandt et al. 2005, 2007

39. All of these are diagonal Nonlocal Continuum Electrostatics: Analytical Solution for Sphere • Solve each mode independently and presto! • Note: This is not about matching interior and exterior expansions--unlike the Kirkwood solution for local model • This decomposition may provide further analytical insights (e.g., eigenvectors of reaction-potential operator) For sphere, these operators share a common eigenbasis: spherical harmonics Bardhan and Brune, to be submitted

40. Local theory needs unrealistically large dielectric constants to match experiment! 3 2 Error in pKa value (RMSD) Measured protein dielectric constants suggest  = 2-5 1 0 5 20 40 60 80 Nonlocal Continuum Electrostatics: Charge Burial and the pKa Problem • Understanding charge burial energetics is important! • For protein folding, misfolding (Alzheimer’s), etc. • For two molecules binding (drug-protein, protein-protein, etc.) • For change in environment (pH, temperature, concentration, etc.) Ion or charged chemical group, alone in water Ion or charged chemical group, buried in protein Demchuk+Wade, 1996

41. Nonlocal Continuum Electrostatics: Charge Burial and the pKa Problem • Nonlocal theory with realistic dielectric constant predicts similar energies as (widely successful) local theories with unrealistic dielectric constants! Bardhan, J. Chem. Phys. (in press)

42. Dense BEM Fast BEM Nonlocal Continuum Electrostatics: Fast Solver is a Must for Accurate Studies • O(N2) memory limitation: big discretization errors • O(N) fast solver: only way to get accurate energies Illustration of surface representations for memory-constrained dense and fast BEM

43. Required accuracy Dense methods used previously could not achieve useful accuracy! Nonlocal Continuum Electrostatics: Fast BIE Solver Performance • Time and memory scale linearly in the number of unknowns • Unoptimized code still allows a laptop to solve 10X larger problems than is possible on a cluster • Preconditioning is vital (use diagonal entries of blocks) Bardhan and Hildebrandt, DAC ‘11

44. Nonlocal Continuum Electrostatics:Fast BIE Solver Enables Tests on Proteins • Observe reduced “electrostatic focusing” • Next step: compare to molecular dynamics Local Model Nonlocal Model Bardhan and Hildebrandt, DAC ‘11

45. Nonlocal Continuum Electrostatics: Adding molecular realism “the right way” Extend the space of models that are supported by good theory Derive fast analytical methods for testing the new theories Test on important open questions Build high-performance solvers for realistic, accurate simulations Applied Math Biophysics Nonlocal model Computer science (HPC)

46. Applied Math Biophysics My research Computer science (HPC) Summary: • Improve understanding of existing models • Develop new models on strong foundations • Stringent tests of new models • Identify critical model weaknesses • Explain previously unresolved phenomena • Leverage HPC expertise by re-using computational primitives • “Think computationally” to gain new insights into model development

47. Acknowledgments • Support: • Wilkinson Fellowship at Argonne National Lab • Partial support from a Rush University New Investigator award • Colleagues: • Ridgway Scott (U. Chicago) • Bob Eisenberg, Dirk Gillespie (Rush) • Mala Radhakrishnan (Wellesley) • Nathan Baker (Pacific Northwest Nat’l Lab)