Boundary-element methods in biological science and engineering

1
Boundary-element methods in biological science
and engineering

• Jaydeep P. Bardhan
• Dept. of Electrical and Computer Engineering
• Northeastern University, Boston MA

2
Four Points For Today

• Cellular and molecular biomedical problems also
  need efficient simulation methods.
• Fast BEM solvers represent an appealing approach
  even at the molecular scale!
• Challenge Persuading community to abandon
  beloved ad hoc fast methods for systematic ones.
• Strategy Systematic methods are more flexible as
  we add new physics and address inverse problems.

3
Point 1 Efficient solvers are needed not only
for macroscopic biomedical problems
10-5 m
100 m
10-1 m
10-9 m
10-2 m
Human body
Organs/ tissues
Human cells
Molecules
Electrocardiography (ECG)
Tumor growth
Lowengrub et al. 09
4
but for microscopic ones as well!
10-5 m
100 m
10-1 m
10-9 m
10-2 m
Human body
Organs/ tissues
Human cells
Molecules
Quantum mechanics
Biomolecule electrostatics and hydrodynamics

• Drug binding
• Protein folding
• Cell physiology
• Molecular design

5
A central molecular-scale modeling problem water.

• Biology uses water to control molecular binding,
  protein folding, etc. Binding example is simple

Protein
6
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

Linearized Poisson- Boltzmann equation
7
A Boundary Integral Method For the Poisson
Biomolecule Problem

1. Boundary conditions handled exactly
2. Point charges are treated exactly
3. Meshing emphasis can be placed directly on the
   interface

8
Fast BEM Solvers are Essential

• Solve Axb approximately 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!

Replace quadratic memory and cubic time
requirements with LINEAR requirements!
9
Application-Specific Challenges

• Continuum-solvent dynamics
• Replace water molecules with dielectric
• Calculate forces at each time step and integrate
• Continuum post-processing of molecular dynamics
• Sample structures from explicit-water MD
• Compute average continuum energy from samples
• Electrostatic component analysis
• Compute each atoms interaction with every other
• Useful in drug design and protein engineering!

10
Communitys Solution Fast Ad Hoc Models

• Up to 100X faster than solving Poisson
• Define effective (nonphysical) parameters
• Plug in to ad hoc (nonphysical) formula

11
Generalized Born theory

• Can give blatantly unphysical results
• exhibits incorrect dependence on dielectric
  constants
• needs all manner of handwaving justifications
  for improvements
• is VERY, VERY popular.

12
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

13
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

14
BIBEE Natural, Rigorous Generalized Born


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
15
Complementary Regimes of Accuracy

• Small molecules reaction potential matrix and
  eigendecomposition (not the integral operator)
• Top right the electric fields induced by several
  eigenvectors of L, at the dielectric boundary
• Charge distributions that generate uniform
  displacement fields are like low-order
  multipoles CFA does well here and P does poorly
• Small eigenvalues are associated with charge
  distributions that generate rapidly varying
  displacement fields these are approximated well
  by P, not CFA

16
Goal Make Fast Models More Rigorous

• Many advantages for chemists/biophysicists
• Enables systematic model improvement
• Prove approximation properties
• Leverage existing fast, scalable algorithms
• Can add better physics as we learn them
• Natural coupling to inverse problems

17
We really want to approximate the dominant modes
of the integral operator.

• The integral operator has to be split into two
  terms
• BIBEE approximates Es 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

18
Mathematical Rigor Enables Systematic
Improvements
Snapshots from MD
BIBEE fluctuations track actual ones very closely
possible applications in uncertainty
quantification
Mean absolute error 4 !

• This effective parameter is expected to be
  rigorously determined by approximating protein as
  ellipsoid (OnufrievSigalov, 06)

BardhanKnepley, J. Chem. Phys. (in press)
19
Goal Make Fast Models More Rigorous

• Many advantages for chemists/biophysicists
• Enables systematic model improvement
• Prove approximation properties
• Leverage existing fast, scalable algorithms
• Can add better physics as we learn them
• Natural coupling to inverse problems

20
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)
21
The Reaction-Potential Matrix

• A weighted combination of charge distributions in
  the solute molecule produces a weighted
  combination of the individual responses
• The canonical basis is the natural, atom-based
  point of view
• We can also use the eigenvector basis for
  analysis!
• In comparing models we dont just have to use the
  total electrostatic solvation free energy

22
Reaction-Potential Operator Eigenvectors Have
Physical Meaning

• Eigenvectors from distinct eigenvalues are
  orthogonal
• Eigenvectors correspond to charge distributions
  that do not interact via solvent polarization
  (this confuses chemists)
• If an approximate method generates a solvation
  matrix , its eigenvectors should line up
  well with the actual eigenvectors, i.e.

i j
23
BIBEE in Separable Geometries

• For half-spaces, spheres, ellipsoids, BIBEE
  exactly reproduces actual eigenvectors.
• Proof for spheres, ellipsoids use appropriate
  harmonics
• Question for future What about near separable
  geometries?

Bardhan and Knepley, 2011
24
BIBEE Is An Accurate, Parameter-Free Model
Snapshots from MD

• Peptide example

Met-enkephalin
BIBEEs stronger diagonal appearance indicates
superior reproduction of the eigenvectors of the
operator.
All models look essentially the same here.
25
Goal Make Fast Models More Rigorous

• Many advantages for chemists/biophysicists
• Enables systematic model improvement
• Prove approximation properties
• Leverage existing fast, scalable algorithms
• Can add better physics as we learn them
• Natural coupling to inverse problems

26
Pre-corrected FFT Algorithm

• Potential calculation is a convolution.
• Convolutions are cheap in frequency space
• Greens function independent! (Laplace,
  Helmholtz, Stokes, etc.)

1. Project charges to grid
2. Point-wise multiplication in frequency space
3. Interpolate grid potentials
4. Pre-correct so that local interactions are
   accurate

Phillips and White (1997)
27
Higher-order Protein BEM
A geometry representative of a protein
The barnase-barstar protein complex
Bardhan Altman et al., 2007 Altman Bardhan,
White, Tidor 2009
28
Develop scalable protein simulations with leaders
in parallel computing FMM

• 760-node GPU cluster Degima

Parallel GPU FMM code
Cost of cluster US 420,000 Sustained 34.6
Tflops Performance/price 80 Mflops/
Application to proteins with PetFMM code
of Yokota, Cruz, Barba, Knepley, Hamada
29
Scalable algorithms enable bigger science
Lysozyme 2K atom charges, 15K surface charges
1000 lysozyme molecules model of a concentrated
protein solution

• How do proteins work in the crowded environment
  of the cell?

Yokota, Bardhan, et al. 2009
30
Goal Make Fast Models More Rigorous

• Many advantages for chemists/biophysicists
• Enables systematic model improvement
• Prove approximation properties
• Leverage existing fast, scalable algorithms
• Can add better physics as we learn them
• Natural coupling to inverse problems

31
We are still adding physics to our models.
Classical modeling one can assume the model is
right!
All simulate same thing!
Circuit simulation Maxwell equations Solid
mechanics elasticity Airplane simulation
Navier-Stokes
Bio-modeling All models are wrong, some are
useful
Diverse set of flawed models. To avoid flaws,
use expert insight. New models are always
evolving! We have to connect multiple models
(uncertainty quantification).
These are just the models associated with the
molecular scale!!
--George Box
32
Adding physics to the continuum model using
nonlocal dielectric theory

• KNOWN weaknesses of Poisson model
• Linear response assumption
• Caveat Nonlinear dielectrics ARE important for
  some molecules!
• Violates continuum-length-scale assumption
• Water molecules have finite size
• Water molecules form semi-structured networks

Test with all-atom molecular dynamics
Relatively small deviation!
yx denotes exactly linear response
Nina, Beglov, Roux 97
33
Nonlocal Continuum ElectrostaticsLorentzian
Model and Promising Tests

• Nonlocal response
• Now
• Integrodifferential Poisson equation

Single parameter fit for ? gives much better
agreement with experiment!!
A. Hildebrandt et al. 2004
34
Nonlocal Continuum ElectrostaticsReformulation
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
Nonlocal Continuum Electrostatics Purely BIE
Formulation

• Three surface variables, two types of Greens
  functions, and a mixed first-second kind problem
• The derivation uses double reciprocity theory,
  which can be applied to nonlinear problems as
  well!
• Have derived exact solution for charges in a
  sphere

Hildebrandt et al. 2005, 2007
36
Just as fast, but now with better physics!

• Unoptimized code still allows a laptop to solve
  10X larger problems than is possible on a cluster
  with dense methods
• Current work comparing to molecular dynamics
  simulations

Bardhan and Hildebrandt, DAC 11
37
Nonlocal Continuum Electrostatics Charge Burial
and the pKa Problem

• Understanding charge burial energetics is
  important!
• For protein folding, misfolding (Alzheimers),
  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
DemchukWade, 1996
38
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)
39
A Common Framework for Multiple Models
BIBEE provides a unifying, scalable approach to
testing and extending new physics.
GB-like fast nonlocal approximate model
Improved GB models
Fast GB-like nonlinear approximations
Explain Coulomb-field approx.
Full nonlinear PB via boundary-integrals
Coupling to fast, scalable algorithms
Advanced PB models (Bikerman, etc.)
Dynamics hybrid explicit/implicit, and fully
implicit
Popular quantum methods couple to exactly our
Poisson problem (polarizable continuum model)
40
Goal Make Fast Models More Rigorous

• Many advantages for chemists/biophysicists
• Enables systematic model improvement
• Prove approximation properties
• Leverage existing fast, scalable algorithms
• Can add better physics as we learn them
• Natural coupling to inverse problems

41
The Value of Systematic Approximations in Inverse
Problems Biomolecule Design

• The electrostatic contribution to binding is
• A total of three simulations is needed.

42
Electrostatic Optimization of BiomoleculesApplic
ations in Analysis and Design

• E. coli chorismate mutase inhibitors
• Analyzed by Kangas and Tidor
• Suggested substitution experimentally verified
  result is the tightest-binding inhibitor yet
  known
• Barnase/barstar protein complex
• Tight-binding complex
• Optimal charge distribution closely matches
  wild-type charge distribution

43
Challenge Optimization is SLOW.
10 min/simulation 2000 simulations
(protein) 2 CPU weeks!!
44
A Novel Method The Reverse-Schur Approach

• For these PDE constraints, we really only need to
  solve multiple systems simultaneously
• The unconstrained problem is therefore
• Constraints can be handled using standard methods
  (Lagrange multipliers, etc.)

45
New Approach is Dozens to Hundreds of Times
Faster, but Formally Exact
Method scales comparably with normal
PDE-constrained approaches
Formally exact calculation
10 min/simulation 20 min/optimization (no
matter how many charges!)
Bardhan et al., 2004 Bardhan et al., 2005
Bardhan et al., 2007 Bardhan et al., 2009
46
BIBEE as Inverse Problem Regularizer

• Regularization can be performed using
  approximate penalty functions
• No linear solve Accurate but 10-20X faster than
  simulation!!
• BIBEE/P captures small eigenvalues very
  accurately ? identify number of directions to
  penalize

47
Application Cyclin-Dependent Kinase 2 and
Inhibitor
PDE-constrained optimization is almost 200 times
faster for this small molecule
Anderson, et al. 2003 (not exactly the optimized
ligand)
Bardhan et al., J. Chem. Theory Comput. (2009)
48
Summary Pushing On All Dimensions
2. Add Realism But Preserve Speed
1. Fast, Scalable Numerical Methods
3. Solve Inverse Problems in Design
4. Unify Theories For New Science
49
Four Points For Today

• Cellular and molecular biomedical problems also
  need efficient simulation methods.
• Fast BEM solvers represent an appealing approach
  even at the molecular scale!
• Challenge Persuading community to abandon
  beloved ad hoc fast methods for systematic ones.
• Strategy Systematic methods are more flexible as
  we add new physics and address inverse problems.

50
Collaborators and Acknowledgments

• Fast methods Michael Altman (Merck), Matt
  Knepley (U. Chicago), Rio Yokota (King Abdullah
  University of Science and Technology), Lorena
  Barba (Boston U.), Tsuyoshi Hamada (Nagasaki U.)
• Nonlocal continuum theory Andreas Hildebrandt
  (Johannes Gutenberg U., Mainz), Peter Brune,
  David Green (SUNY Stony Brook)
• Fast optimization Michael Altman, Bruce Tidor
  (MIT), Jacob White (MIT), Jung Hoon Lee (Merck),
  Sven Leyffer (Argonne) , Steve Benson (Argonne),
  David Green, Mala Radhakrishnan (Wellesley)
• Approximation method Matt Knepley, Mihai
  Anitescu (Argonne), Mala Radhakrishnan
• Support from
• Department of Energy (DOE) Computational Science
  Graduate Fellowship (CSGF)
• Wilkinson Fellowship in Math and Computer Science
  Division of Argonne National Lab
• NIH Technology Development (EUREKA)
• Rush New Investigator Award