Computational Methods for Biomolecular Diffusion and Electrostatics

1
Computational Methods for Biomolecular Diffusion
and Electrostatics

• Nathan A. Baker
• Department of Biochemistry and Molecular
  Biophysics
• Center for Computational Biology
• Washington University in St. Louis
• BME 140
• November 10, 2003

2
Overview

• Computational electrostatics
• Background
• Applications
• Advances
• Computational diffusion
• Background
• Applications
• Advances
• Outlook

3
Introduction to biomolecular electrostatics

• Highly relevant to biological function
• Important tools in interpretation of structure
  and function
• Electrostatics pose one of the most challenging
  aspects of biomolecular simulation
• Long range
• Divergent
• Existing methods limit size of systems to be
  studied

Acetylcholinesterase
Fasciculin-2
4
Implicit solvent simulations background

• Solute typically only accounts for 5-10 of atoms
  in explicit solvent simulation
• Implicit methods
• Solvent treated as continuum of infinitesimal
  dipoles
• Ions treated as continuum of charge
• Some deficiencies
• Polarization response is linear and local
• Mean field ion distribution ignores fluctuations
  and correlations
• Apolar effects treated by various, heuristic
  methods

5
Modeling biomolecule-solvent interactions

• Solvent models
• Explicit
• Molecular dynamics
• Monte Carlo
• Integral equation
• RISM
• 3D methods
• DFT
• Primitive
• Poisson equation
• Phenomenological
• Generalized Born
• Modified Coulombs law

• Ion models
• Explicit
• Molecular dynamics
• Monte Carlo
• Integral equation
• RISM
• 3D methods
• DFT
• Field theoretic
• Poisson-Boltzmann
• Extended PB, etc.
• Phenomenological
• Generalized Born
• Debye-Hückel

Level of detail
Computational cost
6
Implicit solvent simulations

• Free energy evaluations
• Usually based on static solute structures or
  small number of conformational snapshots
• Solvent effects included in
• Implicit solvent electrostatics
• Surface area-dependent apolar terms
• Useful for
• Solvation energies
• Binding energies
• Mutagenesis studies
• pKa calculations

Animation courtesy of Dave Sept
7
Implicit solvent simulations

• Stochastic dynamics
• Usually based on Langevin or Brownian equations
  of motion
• Solvent effects included in
• Implicit solvent electrostatics forces
• Hydrodynamics
• Random solvent forces
• Useful for
• Bimolecular rate constants
• Conformational sampling
• Dynamical properties

Animation courtesy of Dave Sept
8
Coulombs law

• Simplest electrostatic description
• Coulombs law solution to the Poisson equation
  for
• Homogeneous dielectric
• Point charges
• Infinite space

9
Coulombs law

• Simplest electrostatic description
• Coulombs law solution to the Poisson equation
  for
• Homogeneous dielectric
• Point charges
• Infinite space
• Potential due to multiple charges obtained by
  superposition

10
Coulombs law

• Simplest electrostatic description
• Coulombs law solution to the Poisson equation
  for
• Homogeneous dielectric
• Point charges
• Infinite space
• Potential due to multiple charges obtained by
  superposition
• Simple energy evaluation

11
Debye-Hückel law

• Similar system to Coulombs law solution to the
  Helmholtz equation
• Added mobile ionic species
• Still have

12
Debye-Hückel law

• Similar system to Coulombs law solution to the
  Helmholtz equation
• Added mobile ionic species
• Still have
• Superposition
• Simple energy evaluation

13
Debye-Hückel law

• Similar system to Coulombs law solution to the
  Helmholtz equation
• Added mobile ionic species
• Still have
• Superposition
• Simple energy evaluation

14
Generalized Born methods

• Modified version of Coulombs law
• Heuristic
• Uses modified distance function with
• Altered decay (often exponential)
• Effective Born radius
• Born radius adjusted to give desired solvation
  energy for each atom
• Radius dependent on location
• Buried atoms have larger radii than exposed atoms
• Potential for multiple charges obtained by
  superposition

15
Poisson-Boltzmann and Poisson equations
16
Poisson-Boltzmann and Poisson equations
17
Poisson-Boltzmann and Poisson equations
18
Poisson-Boltzmann and Poisson equations
19
Poisson-Boltzmann and Poisson equations
20
The Poisson-Boltzmann equation current solution
methods

• Complicated geometries require numerical
  solutions
• Numerical methods
• Local vs. global basis functions
• Discretization
• Finite domain (usually) with appropriate boundary
  conditions
• PB methods usually use local basis functions
  spatial discretization
• Beware numerical artifacts!
• Convergence of the method
• Inappropriate spacings

21
The Poisson-Boltzmann equation current solution
methods
22
The Poisson-Boltzmann equation current solution
methods
23
Free energy evaluation

• Integral of electrostatic potential over solution
  domain
• Includes contributions from

24
Free energy evaluation

• Integral of electrostatic potential over solution
  domain
• Includes contributions from

25
Free energy evaluation

• Integral of electrostatic potential over solution
  domain
• Includes contributions from

26
Free energy evaluation

• Integral of electrostatic potential over solution
  domain
• Includes contributions from
• Can also be linearized

27
A continuum descriptionof ion solvation

• Born ion model
• Non-polarizable ion
• Point charge
• Higher polarizability medium
• Reaction field effects
• Non-Coulombic potential inside ion due to
  polarization of solvent
• Solvation energy
• Simple model with analytical solutions

28
A continuum descriptionof ion solvation
29
A continuum descriptionof ion desolvation

• Two Born ions at varying separations
• Solve Poisson equation at each separation
• Increase in energy as water is squeezed out of
  interface
• Desolvation effect
• Less volume of polarized water
• Important points
• Non-superposition of Born ion potentials
• Reaction field causes repulsion at short
  distances
• Dielectric medium focuses field

30
A continuum descriptionof ion desolvation
31
Non-specific salt effects screening

• Lots of types of non-specific ion screening
• Variable solvation effects (Hofmeister)
• Ion clouds damping electrostatc potential
• Changes in co-ion and ligand activity
  coefficients
• Condensation
• Not all ion effects are non-specific!
• Generally reduces effective range of
  electrostatic potential
• Shown here for acetylcholinesterase
• Illustrated by potential isocontours
• Observed experimentally in reduced binding rate
  constants

32
Non-specific salt effects screening
33
Electrostatic influenceson ligand binding

• Examine inhibitor binding to protein kinase A
• Part of drug design project by McCammon and
  co-workers
• Illustrates how electrostatics governs
  specificity and affinity
• Look at complementarity between ligand and
  protein electrostatics
• Verify with experimental data (relative binding
  affinities)
• Use to guide design of improved inhibitors

34
Electrostatic influenceson ligand binding
35
Electrostatic influenceson ligand binding
36
Poisson-Boltzmann equationnew solution methods

• Parallel multilevel adaptive finite element
  techniques
• Adaptivity places computational effort where
  needed
• Trivially parallel
• Additional overhead due to unstructured mesh
• Currently suited to smaller systems or small
  molecule binding
• Parallel focusing methods
• Very fast multigrid-based solver
• Trivially parallel
• Works at all scales
• No adaptivity

37
Parallel focusing algorithm

• Given the problem data and P processors of a
  parallel machine
• Each processor i 1, , P
• Obtains a coarse solution over the global domain
• Subdivides the global domain into P subdomains,
  each of which is assigned a processor
• Assigns boundary conditions to a fine
  discretization of its subdomain using the coarse
  global solution
• Solves the equation on its subdomain
• A master processor collects observable data from
  other processors and controls I/O

38
Parallel focusingbenefits

• Loosely coupled focusing calculations
• Trivially parallel algorithm
• No load balancing issues
• Simple implementation
• Leverage existing, highly optimized multigrid
  solvers
• Simple force and observable evaluation

39
Parallel focusing application to ribosomes

• Ribosome central to protein synthesis machinery
• Target for several pharmaceuticals
• Nucleoprotein composition make it computationally
  challenging
• Composed of two subunits (large and small)
• 30S consists of 88,000 atoms and roughly 200 Å
  cube
• 50S consists of more than 95,000 atoms and
  roughly 200 Å cube
• Function involves several interesting features
• Protein-nucleic acid association
• Protein-protein association
• Conformational changes
• Salt dependence (type and quantity)
• Solved on 343 processors of Blue Horizon to 0.41
  Å (30S) and 0.43 Å (50S) resolution

40
Parallel focusing application to ribosomes
41
Parallel focusing application to ribosomes
42
Parallel focusing ribosome-antibiotic
interactions

• Determine binding energies between 30S ribosomal
  subunit and aminoglycoside antibiotics

43
Parallel focusing ribosome-antibiotic
interactions

• Determine binding energies between 30S ribosomal
  subunit and aminoglycoside antibiotics

44
Parallel focusing ribosome-antibiotic
interactions

• Determine binding energies between 30S ribosomal
  subunit and aminoglycoside antibiotics
• Excellent fit to data 0.78 0.13 slope with
  small antibiotics, 0.95 0.19 slope without

45
Parallel focusing ribosome-antibiotic
interactions

• Determine binding energies between 30S ribosomal
  subunit and aminoglycoside antibiotics
• Excellent fit to data 0.78 0.13 slope with
  small antibiotics, 0.95 0.19 slope without
• Suggests importance of basic functional groups on
  Ring IV

46
Parallel focusingapplication to microtubules

• Important cytoskeletal components structure,
  transport, motility, division
• Typically 250-300 Å in diameter and up to
  millimeters in length
• Computationally difficult due to size (1,500
  atoms/Å ) and charge (-4.5 e/Å)
• Solved LPBE at 150 mM ionic strength on 686
  processors for 600 Å-long, 1.2-million-atom
  microtubule
• Resolution to 0.54 Å for largest calculation
  quantitative accuracy

47
Parallel focusingapplication to microtubules
48
Parallel focusingapplication to microtubules
49
Parallel focusingmicrotubule stability and
assembly

• Performed series of calculations on tubulin
  dimers and protofilament pairs
• Poisson-Boltzmann electrostatics and SASA apolar
  energies
• Observed 7 kcal/mol stronger interactions between
  protofilaments than within
• Determined energetics for helix properties
  predict correct minimum for experimentally-observe
  d A (52 Å) and B (8-9 Å) lattices

50
Parallel focusingmicrotubule stability and
assembly
51
Parallel focusingmicrotubule stability and
assembly
52
Quantitative analyses of electrostatic potentials

• Current progress
• TAQT
• Under development
• Computes topological properties of scalar data
  (Betti numbers, Morse complexes)
• Computes contour spectrum of data
• Similarity comparisons
• Various norms available
• Datasets can be compared as is or after
  multiresolution analysis (multipole expansion,
  downsampling)
• Under testing to be released as part of APBS
• Future work
• Integrate with electrostatics pipeline
• Large-scale comparisons within protein functional
  and structural classes
• Cache data for querying and/or comparisons with
  new potentials

53
The fine print (Part 1)Breakdown of implicit
solvent models

• Is the water near a zwitterionic bilayer a
  featureless dielectric continuum?
• Examine behavior of TIP3P water around POPC
  bilayer using molecular dynamics
• Systematically replace each layer of explicit
  solvent with dielectric continuum and calculate
  potential
• Average results over MD trajectory
• Determine
• Nature of water at membrane surface
• Discrepancies between implicit and explicit
  electrostatics

54
Water around membrane systems results

• 4 layers of water significantly different than
  bulk
• Dielectric response (orientational fluctuations)
  of water near membrane much lower
• Membrane potential dramatically affected by water
  polarization
• Conclusion Water near membrane has
  significantly different structure and response
  properties from bulk
• Future work non-neutral membranes and mobile
  ions

55
The fine print (Part 2)Breakdown of the
implicit ion model

• Canonical electrostatics test case
  non-polarizable hard sphere with point charge
• Mild test of continuum electrostatics (no
  condensation, etc.)
• Good agreement with mean-field results at low
  mobile ion concentrations
• Agreement deteriorates with
• Increasing macroion charge
• Increasing mobile ion concentration

56
The fine print (Part 2)Breakdown of the
implicit ion model
57
Ions around DNA

• NPBE simulations predict 2000 M Mg near DNA
  surface for 20 mM MgCl2, 150 mM NaCl solution
• Ran GCMC calculations for similar conditions
  around 20 bp D-DNA
• Observe max density of 10 M Mg in GCMC
  simulations
• Much lower than PBE results
• Implicit solvent model questionable
• Reproduce expected condensation behavior
• 88 charge compensation at 17 Å radius
• Divalent cations affect charge screening

58
Ions around DNA
59
Ions around DNA
60
A plug APBS software

• APBS Adaptive Poisson-Boltzmann Solver
• Runs on machines from Macs to Crays
• Currently available for download from
  http//agave.wustl.edu/apbs
• Further reading many of todays examples
  available as APBS scripts

61
ElectrostaticsSummary and outlook

• Poisson-Boltzmann equation an implicit solvent
  method for biomolecular electrostatics

62
ElectrostaticsSummary and outlook

• Poisson-Boltzmann equation an implicit solvent
  method for biomolecular electrostatics
• New algorithms enable fast and scalable solution
  of the PBE

63
ElectrostaticsSummary and outlook

• Poisson-Boltzmann equation an implicit solvent
  method for biomolecular electrostatics
• New algorithms enable fast and scalable solution
  of the PBE
• Applied to large biomolecular systems

64
ElectrostaticsSummary and outlook

• Poisson-Boltzmann equation an implicit solvent
  method for biomolecular electrostatics
• New algorithms enable fast and scalable solution
  of the PBE
• Applied to large biomolecular systems
• Integrate into pipeline for automated
  electrostatic calculations and characterization

65
ElectrostaticsSummary and outlook

• Poisson-Boltzmann equation an implicit solvent
  method for biomolecular electrostatics
• New algorithms enable fast and scalable solution
  of the PBE
• Applied to large biomolecular systems
• Integrate into pipeline for automated
  electrostatic calculations and characterization
• Where do we go from here?
• Use electrostatic pipeline for biomolecular
  comparison
• Investigate limitations of PBE framework
• Integrate into multiscale methods

66
Introduction to biomolecular diffusion

• Highly relevant to biological function
• Passive transport
• Biomolecular association
• Catalysis
• Multiple levels of detail
• Solvent
• Conformational changes
• Chemistry
• Electrostatics play an important role!

67
Methods for biomolecular diffusion

• Discrete methods
• Usually employed reduced models
• Integrate stochastic equations of motion
• Pros offer explicit atomic detail and
  incorporation of other stochastic phenomena
• Cons slow convergence, difficult to scale

68
Methods for biomolecular diffusion

• Discrete methods
• Usually employed reduced models
• Integrate stochastic equations of motion
• Pros offer explicit atomic detail and
  incorporation of other stochastic phenomena
• Cons slow convergence, difficult to scale
• Continuum methods
• Reduced models solvent and ligands
• Deterministic
• Pros good convergence, highly scalable, provide
  connections to other continuum mechanics behavior
• Cons lack of detail/correlation in diffusing
  species

69
Smoluchowski equation
Non-reactive boundary
Reactive boundary
Outer boundary
70
Smoluchowski equation
71
Smoluchowski equation
72
Smoluchowski equation
Non-reactive boundary
Reactive boundary
Outer boundary
73
Smoluchowski equation
74
Observable reaction rates

• Measure flux through reactive surface(s)
• Directly related to reaction rate observed via
  stochastic methods (cf. Zhou HX, J Chem Phys 92
  (5) 3092-3095, 1990)
• As with BD, measure rate as function of
  mutation, environmental conditions, diffusing
  species, etc.

75
Our model

• Boundaries
• Domain inner boundaries determined by
  biomolecular surface
• Reactive portion determined by experimental
  knowledge of active/binding sites
• Outer boundary set at large distance bulk
  density set to unit concentration

76
Our model

• Boundaries
• Domain inner boundaries determined by
  biomolecular surface
• Reactive portion determined by experimental
  knowledge of active/binding sites
• Outer boundary set at large distance bulk
  density set to unit concentration
• Potential of mean force (W)
• Electrostatic forces determined from PBE
• Currently not coupled to diffusion

77
Our model

• Boundaries
• Domain inner boundaries determined by
  biomolecular surface
• Reactive portion determined by experimental
  knowledge of active/binding sites
• Outer boundary set at large distance bulk
  density set to unit concentration
• Potential of mean force (W)
• Electrostatic forces determined from PBE
• Currently not coupled to diffusion
• Kinetic variables
• Diffusion constant set from Stokes-Einstein
• Steady-state and time-dependent versions solved

78
Solution techniques

• Discretization with adaptive finite element basis
• Steady-state and per-time-step solutions
• Iterative asymmetric solvers
• Multilevel methods
• Time-dependent solutions
• Fixed finite element mesh
• Variety of ODE solvers (backward Euler)
• Electrostatic potential
• APBS PBE solutions
• Multiple grids/multiple resolutions

79
Solution techniques finite element methods

• Discretization on adapted finite element mesh

80
Solution techniques finite element methods

• Discretization on adapted finite element mesh

81
Solution techniques finite element methods

• Discretization on adapted finite element mesh
• Linear systems solved by multilevel method using
  algebraic hierarchy of operators
• Ak1 PTk Ak Pk
• Nonlinear systems solved with inexact-damped
  Newton

82
Solution techniques finite element methods

• Discretization on adapted finite element mesh
• Linear systems solved by multilevel method using
  algebraic hierarchy of operators
• Ak1 PTk Ak Pk
• Nonlinear systems solved with inexact-damped
  Newton
• Parallel solution via Bank-Holst method
• Optimal to near-optimal numerical performance

83
Solution techniquesmesh generation

• The devils in the details
• Adaptive mesh generation
• C. Bajaj group (UT Austin)
• Dual Contouring Method
• Mesh quality refinement
• Adaptive mesh refinement
• Error estimation
• Geometry-based refinement
• Simplex subdivision
• Result nested meshes

84
Solution techniquesmesh generation

• The devils in the details
• Adaptive mesh generation
• C. Bajaj group (UT Austin)
• Dual Contouring Method
• Mesh quality refinement
• Adaptive mesh refinement
• Error estimation
• Geometry-based refinement
• Simplex subdivision
• Result nested meshes

85
Validationuniformly-reactive sphere

• Assuming a spherical cow
• Good test case with analytical results for
  comparison
• System
• Spherical fixed molecule
• Coulombic potential
• Water friction coefficient
• Results
• Good agreement with analytical results
• Fast 3-5 minutes per run
• Unstable (as with analytical case) at low
  dielectric screening

86
Validationuniformly-reactive sphere
87
Acetylcholinesterase kinetics

• Using mAChE structure prepared by McCammon group
• Binds variety of positively-charged ligands at
  diffusion-controlled limit (107 to 1010 M-1 s-1)
• Extensive kinetic measurements
• Experimental
• Brownian dynamics
• Good test case for validation in realistic
  geometries

88
mAChE mesh generation

• Several difficult issues
• Typical bumpy molecular surface
• Highly charged
• Long, narrow gorge
• Span very large domain adaptively
• Final meshes extremely high quality

89
mAChE mesh generation

• Several difficult issues
• Typical bumpy molecular surface
• Highly charged
• Long, narrow gorge
• Span very large domain adaptively
• Final meshes extremely high quality

90
mAChE mesh generation

• Several difficult issues
• Typical bumpy molecular surface
• Highly charged
• Long, narrow gorge
• Span very large domain adaptively
• Final meshes extremely high quality
• High simplex size ratio between inner and outer
  boundaries
• No error-based refinement (yet)

91
mAChE kinetics

• Runs performed to compare with McCammon group BD
  results and experimental data
• Calibrate reactive surface on experimental data
• Determine surface reaction criteria
• Ligand modeled as 2 Å sphere (TMA) neither ATCh
  nor TMTFA
• BD results log(kon) 11 at 670 mM ionic
  strength
• Interesting changes in reaction rate near ASP74
• Like BD, probably currently most useful for
  relative rate changes
• however, much faster 10 min compared to 60-90
  min

92
Large-scale applicationsneuromuscular junctions

• Collaboration with J. A. McCammon, M. Holst, M.
  Ellisman (UCSD)
• ACh diffusion in a neuromuscular juction
• Release from vesicle
• Hydrolysis by AChE
• Binding to AChR
• Little or no molecular details only gross
  morphology

93
Large-scale applicationsneuromuscular junctions

• Collaboration with J. A. McCammon, M. Holst, M.
  Ellisman (UCSD)
• ACh diffusion in a neuromuscular juction
• Release from vesicle
• Hydrolysis by AChE
• Binding to AChR
• Little or no molecular details only gross
  morphology
• Qualitative agreement with mouse MEPCs for fast-
  and slow-twitch muscles

94
Summary and outlook

• Develop robust continuum diffusion solver and
  software
• Improved mesh generation
• Faster solvers
• Publicly-available code
• Integrate other continuum mechanics behavior
• Electrostatics (PNP)
• Elasticity
• Fluid dynamics
• Multiscale modeling

95
Summary and outlookmultiscale modeling

• Examine the effect of
• Molecular information on large-scale systems
• Large-scale changes on molecular structure and
  function
• Applications
• Synaptic transmission
• Biomolecular elasticity

96
Summary and outlookmultiscale modeling

• Examine the effect of
• Molecular information on large-scale systems
• Large-scale changes on molecular structure and
  function
• Applications
• Synaptic transmission
• Biomolecular elasticity

97
Acknowledgements

• APBS
• Andy McCammon
• Mike Holst
• Microtubules
• Dave Sept
• Ribosome
• Chiansan Ma
• Simpson Joseph
• Membrane
• Jung-Hsin Lin, Yuhua Song
• Continuum diffusion
• Yuhua Song, Tongye Shen, Andy McCammon, Jessica
  Zhang, Chandrajit Bajaj

Electrostatic informatics Chandrajit Bajaj, Jens
Nielsen, Todd Dolinsky, Jerry Greenberg, Zaiqing
Xu Computational resources NPACI/SDSC/UT, NBCR,
Keck Centers, PSC, IBM Funding NPACI, Wash U
Animation by Jerry Greenberg