Transcript How Well Does Poisson-Boltzmann Implicit Solvent Compare

Poisson-Boltzmann Molecular Dynamics:
Theory and Algorithms
Ray Luo
Molecular Biology and Biochemistry
University of California, Irvine
1
Different levels of abstraction:
Approximations of biomolecules
•
Quantum description: electronic & covalent structure
•
Atom-based description: non-covalent interactions
•
Residue-based/coarse-grained description: overall
motion/properties of a biomolecule
2
Intermolecular forces
Intermolecular Forces, A.J. Stone
3
Biomolecules on computer:
Classical molecular mechanics
Kb (b  b0 )2
K (  0 )2
Bonded
Potential
Energy
K 
2
Electrostatic
Nonbonded
V 1 cos(n   )
0
n
n
qi q j
 rij
Repulsion-dispersion
A
B

rij12
rij6
4
Challenges in biomolecular simulations:
Atomistic representation
•
•
Realistic water environment
Long-range interactions
•
•
Periodic boundary
How to avoid O(n2)?
5
Challenges in biomolecular simulations:
Time scales are in the 109 time steps
Salt bridge population
0.5
0.4
0.3
0.2
0.1
1
2
3
4
5
6
7
8
9
Times (x10ns)
Multiple trajectories, often as many as 10s to 100s, are needed
6
Explicit solvent and implicit solvent:
Removing solvent degrees of freedom
P(ru , rv ) 
exp[U (ru , rv )]
 drudrv exp[U (ru , rv )]
exp[W (ru )]   drv exp[ U (ru , rv )]
P (ru ) 
exp[  W (ru )]
 dru exp[ W (ru )]
ru: solute coordinates; rv: solvent coordinates
7
Continuum solvation approximations
• Homogenous structureless solvent distribution
• Solute geometry (shape/size) influence in solvent
density is weak in solvation free energy calculation
• Solvation free energy can be decomposed into
different components
W  Wpol  Wnpol
Wnpol  Wrep  Watt
8
Polar solvation
   (r) (r)   (r)   ni0qi exp[ qi (r)]
Dielectric
constant
Charge density
p+
Charge of salt ion in solution
Electrostatic potential
+
+
+
s
9
Nonpolar solvation
Wnpol  Wrep  Watt
Wrep : Estimated with surface (SES/SAS) or volume (SEV/SAV)
Wrep   A  c
Wrep  pV  c
uv
Watt: Approximated by Uatt
(D. Chandler and R. Levy)
Nu
uv
Uatt
   neff Vatt (r)d 3r
i 1
10
Is Continuum Approximation Sufficient?
I. Polar Solvation
11
Explicit solvent (TI)
•
TIP3P water model. Periodical
Boundary Condition. Particle Mesh
Ewald, real space cutoff 9Å.
•
NPT ensemble, 300K, 1bar. Preequilibrium runs at least 4 ns and until
running potential energy shows no
systematic drift.
•
All atoms restrained to compare with
PB calculations on static structures
•
25 λ’s with simulation length doubled
until free energies change less than
0.25kcal/mol (up to 320ps
equilibration/production per λ needed).
•
Thermodynamic Integration:
G  
1
0
dH( )
d
d

12
Implicit solvent (PB)
• Final grid spacing 0.25 Å. Two-level focusing
was used. Convergence to 10-4.
• Solvent excluded surface. Harmonic dielectric
smoothing was applied at dielectric boundary.
• Charging free energies were computed with
induced surface charges.
• (110+110 snapshots) × 27 random grid origins
were used.
• Cavity radii were refitted before comparison
Linearized Poisson-Boltzmann Equation:
ε= 80
where
13
Fitting quality: Polar solvation free energies
Correlation Coefficient:
0
Root Mean Square Deviation:
0.33 kcal/mol
AMBER/TIP3P Error (wrt Expt):
1.06 kcal/mol
AMBER/PB Error (wrt Expt):
0.97 kcal/mol
(neutral side chain analogs)
Gelec by PB (kcal/mol)
0.99995
-20
-40
-60
-80
-100
-100
-80
-60
-40
-20
0
Gelec by TI (kcal/mol)
Tan et al, JPC-B, 110, 18680-18687, 2006
14
Salt-bridge charging free energies
(a) Tested salt bridge with atom ids.
(b) PEPenh, a 16mer helix from1enh.
(c) ENH, (1enh, ~50 aa).
(d) P53a, (1tsr, ~200 aa)
ARG154-GLU76 on p53.
(a) P53b, ARG178-GLU190 on p53.
Tan and Luo, In Prep.
15
Salt-bridge charging free energies
Gele_exp(kcal/mol)
-20
-40
-60
-80
p53
-100
PEP
enh
PEP
ENH
P53
-120
-120 -100 -80 -60 -40
Gele_imp(kcal/mol)
-20
Tan and Luo, In Prep
16
Is Continuum Approximation Sufficient?
II. Nonpolar Solvation
17
Explicit solvent (TI)
• TIP3P water model. Periodical Boundary Condition. Particle
Mesh Ewald, real space cutoff 9Å.
• NPT ensemble, 300K, 1bar. Pre-equilibrium runs with neutral
molecules for at least 8 ns and until running potential energy
shows no systematic drift.
• All atoms restrained to compare with single-snapshot
calculations in implicit solvent.
• Thermodynamic Integration:
G  
1
0
dH( )
d
d

• 60 λ’s with simulation length doubled until free energies
change less than 0.25kcal/mol (160ps equilibration or
production per λ needed).
Tan et al, JPC-B, 111, 12263-12274, 2007
18
(A) SES
CC: 0.997
RMSD: 0.30kcal/mol
RMS Rel Dev: 0.026
(B) SEV
CC: 0.985.
RMSD: 0.69kcal/mol
RMS Rel Dev: 0.082
(C) SAS
CC: 0.997
RMSD: 0.30kcal/mol
RMS Rel Dev: 0.026
(D) SAV
CC: 0.998.
RMSD: 0.27kcal/mol
RMS Rel Dev: 0.022
Grep_exp (kcal/mol)
Fitting Quality:
Nonpolar repulsive free energies
25
A:SES
B:SEV
20
15
10
5
0
D:SAV
C:SAS
20
15
10
5
0
0 5 10 15 20 0 5 10 15 20 25
Grep_imp (kcal/mol)
Tan et al, JPC-B, 111, 12263-12274, 2007
19
Fitting quality:
Nonpolar attractive free energies
CC: 0.9995
RMSD: 0.16kcal/mol
RMS Rel Dev: 0.01
Gatt_exp (kcal/mol)
0
-5
-10
-15
-20
-25
-25
Tan et al, JPC-B, 111, 12263-12274, 2007
-20 -15 -10 -5
Gatt_imp (kcal/mol)
Error bars too small to be seen
0
20
Nonpolar solvation free energies of TYR
(a) Tested side chain with atom ids.
(b) PEPα, a 17mer helix from 1pgb.
(c) PEPβ, a 16mer hairpin from 1pgb.
(d) PGB, 1pgb, ~50 aa.
(e) P53, 1tsr, ~200 aa.
Tan and Luo, In Prep.
21
Nonpolar attractive free energies
-2
CC: 0.983
RMSD: 0.29 kcal/mol
RMS Rel Dev: 0.035
Gatt_exp (kcal/mol)
-4
-6
-8

PEP

PEP
PGB
P53
-10
-12
-12
-10
-8
-6
-4
Gatt_imp (kcal/mol)
-2
Tan and Luo, In Prep.
Error bars too small to be seen
22
Nonpolar repulsive free energies
16
12
8
(B) SAV
CC: 0.984
RMSD: 0.53kcal/mol
RMS Rel Dev: 0.053
Tan and Luo, In Prep.
4
Grep_exp (kcal/mol)
(A) SAS
CC: 0.975
RMSD: 2.42kcal/mol.
RMS Rel Dev: 0.55

PEP

PEP
PGB
P53
0
A
-4
12
8
4
0
B
-4
-4
0
4
8
12
Grep_imp (kcal/mol)
16
23
Behaviors of Two Estimators for
TYR Side-Chain Conformations
SAS
SAV
0.06
0.3
3
p (kcal/mol/Å )
0.05
2
 (kcal/mol/Å )
0.2
0.1
0.0

PEP

PEP
PGB
P53
AVG_HLX
-0.1
-0.2
0.04

0.03
0.02
0
2
4
6
8 10
Conformations
PEP

PEP
PGB
P53
AVG_HLX
12
0
2
4
6
8 10
Conformations
12
Tan and Luo, In Prep.
24
Continuum solvation approximation
•
Conformation dependent energetics is consistent
between implicit and explicit solvents.
•
Both polar and nonpolar attractive component correlate
very well with TI from short peptides up to proteins of
typical sizes.
•
Repulsive nonpolar component works well from tested
peptides to proteins if the volume estimator is used.
25
Going beyond Fixed Charge Models with
Continuum Electronic Polarization
26
How to include polarization
in implicit solvents?
•
Explicit treatment
Maple, Cao, et al., J Chem Theo Comp, 1:694, 2005.
Schnieders, Baker, et al., J Chem Phys, 126:124114, 2007.
•
Implicit treatment
27
Continuum polarizable force field
•
Relation between P and E
•
Relation between  and ε
P
Solute dielectric constant ε is optimized
•
P is defined within the molecular volume (solvent
excluded volume).
28
Continuum polarizable force filed
Tan and Luo, J Chem Phys, 126:094103, 2007.
Tan, Wang, and Luo, J Phys Chem, 112:7675. 2008.
29
Continuum polarizable force field
•
Advantage: gives us an efficient and self-consistent
approach in treating polar interactions in biomolecular
simulations more satisfactory than existing additive force
fields with implicit solvents.
•
Limitation: lack of atomic-detailed polarization within a
molecular environment. This may be overcome by use of
functional-group-specific dielectric constants.
30
Charge derivation procedure: RESP
Yes
Convergence
No
Tan and Luo, J Chem Phys, 126:094103, 2007.
31
Quantum mechanical field
•
Computation of quantum mechanically
electrostatic field:
1) Optimization with HF/6-31G*
2) Single point with B3LYP/cc-pVTZ
•
PCM was used for modeling polarization
responses to different environments.
32
monomers
0.9
largest error
rmsd
rmsd and largest error (D)
rmsd and largest error (D)
Quality of fit: dielectric constant
0.6
0.3
0.0
2
3
4
5
solute dielectric constant
dimers
largest error
rmsd
1.5
1.2
0.9
0.6
2
3
4
5
solute dielectric constant
Left: 12 monomers in three environments (vacuum, ε = 4, water)
Right: 4 dimers in three environments
atomic radii: UA0
probe radius:1.385Å
33
Fitting statistics for monomers
in vacuo
8
10
solvent  = 4.0
8
solvent  = 78.4
6
4
2
B3LYP/cc-pVTZ (D)
B3LYP/cc-pVTZ (D)
B3LYP/cc-pVTZ (D)
8
6
4
2
polarizable
nonpolarizable
0
polarizable
nonpolarizable
0
0
2
4
6
8
0
2
4
6
8
6
4
2
polarizable
nonpolarizable
0
0
2
6
8
10
MM (D)
MM (D)
MM (D)
4
in vacuo
ε = 4.0
ε = 78.4
rmsd
0.1425
0.2106
0.1759
uavg
0.1176
0.1848
0.1345
correlation
0.9979
0.9995
0.9996
Dipole moments of monomer with charges fitted simultaneously in three environments
34
Unit: Debye
Transferability among conformations
B3LYP/cc-pVTZ (D)
9
6
3
in vacuo
in =4
in =78.4
0
0
3
6
9
MM (D)
rmsd: 0.2799 uavg: 0.2413 correlation: 0.9922
charges fitted simultaneously for both alphaL and c7eq
in three environments
35
Continuum electronic polarization
•
Electronic polarization with a continuum dipole moment
density. The uniform solute dielectric constant is the only
parameter.
•
Performance comparable to ff02 explicit polarizable force
field for tested dipole moments in vacuum.
•
A single set of charges can be used in different
environments and different conformations. The model
transfers well from monomers to dimers.
36
Poisson-Boltzmann Molecular Dynamics
37
Singular Charges in PBE
•
 function in the PBE
  (r) (r)  4 (r)  4  ni0qi exp[ qi (r)]
 (r)   qi  r  ri 
i
•
Challenges
- Large error in potential near singular charges
- Large error in dielectric boundary force
- Self energy between redistributed charges
38
Removal of Charge Singularity
•
Solve the Laplace’s equation for reaction field potential inside and
simultaneously solve Poisson-Boltzmann equation for total potential
outside.
 

•
   RF
Reaction potential is the difference between the total potential
RF    C
•
Coulombic potential, which is defined as
 2C  40
Cai, Q. et al. Journal of Chemical Physics. 2009, 130, 145101.
39
Removal of Charge Singularity
2 


  0
  2 






f

 0


inside
outside
      C










 
 C
n
n
n
On the dielectric boundary
Cai, Q. et al. Journal of Chemical Physics. 2009, 130, 145101.
40
Discontinuous Interface
  (r) (r)  4 (r)  4  ni0qi exp[ qi (r)]
•
Boundary conditions on the discontinuous interface of the
PBE (uniform potential)
- The potential is continuous on the interface
  
- Integrating the PBE and then using the Gauss’s law
give the flux condition


 
 


n
n
41
Harmonic Average (HA)
•
This method enforces the flux conditions in the three
orthogonal directions on the physical interface, i.e.,








x
x
•

 
 


y
y









z
z
The dielectric constant between two grid points that are
in two different regions is a harmonic average of the
two dielectric constants of the two regions.
h
 1

  i  , j, k  
 2
 a  b




Davis and McCammon, Journal of Computational Chemistry. 1991, 12, 909.
42
Immersed Interface Method (IIM)
•
•
A more accurate method for interface treatment for FDM
IIM proposes new equations involving 27 points instead of the
original 7-point finite-difference equations at the points close to the
interface.
27

m 1
•
 (i  im , j  jm , k  km )  f (i, j, k )  C (i, j, k )
m
IIM tries to minimize the local truncation error with the help of
interface conditions.
27
T (i, j, k )    m (i  im , j  jm , k  km )  f (i, j, k )  C (i, j, k )
m 1
LeVeque and Li. SIAM Journal Numerical Analysis. 1994, 31, 1019.
43
IIM + Removal of Singularity
Tested in the Poisson equation:
single particle system, dielectric boundary force
IIM−Singularity
HA−Singularity
HA+Singularity
Max Error
Max Error
Max Error
4
0.000005
0.000287
0.007136
0.25
8
0.000002
0.000137
0.003402
0.25
16
0.000000
0.000064
0.000845
1.00
4
0.004458
0.004006
1.572224
1.00
8
0.000940
0.001264
0.019043
1.00
16
0.000223
0.000376
0.009309
1.50
4
0.916579
1.148297
59.290638
1.50
8
0.100221
0.096707
6.433441
1.50
16
0.009313
0.010432
0.075564
d
1/h
0.25
Wang, J. et al. Chemical Physics Letters. 2009, 468, 112.
44
Dielectric boundary force: Theory

fbnd
 Pi  e  Po  e
1


2
(

E
E


E
)
e
e
(

E
E
)
e
e
(

E
E
)
e
e
i i i
 
i i i
 
 i i i 2 i i  


 e
1 
1
 0
Pi  e 
( i Ei Ei )e e
( i Ei Ei   i Ei2 )e e
( i Ei Ei )e e

4 
2

 0
1
2

( i Ei Ei )e e
( i Ei Ei )e e
( i Ei Ei   i Ei )e e 
2


1
1
1
1

( i Ei2   i Ei2 )e 
( i Ei Ei )e 
( i Ei Ei )e
4
2
4
4
1
1
1
1
Po  e 
( o Eo2   o Eo2 )e 
( o Eo Eo )e 
( o Eo Eo )e
4
2
4
4
 i Ei   o Eo

 Ei  Eo
 E E
o
 i

f bnd

1
4
1
1

2
2
2
2 
(

E


E
)

(

E


E
)
i
i

i
i
o
o

o
o

2
2

45
Dielectric boundary force: Theory

f bnd
  Pi (0)  e  Pi (  )  e    Po (0)  e  Po (  )  e   0

f bnd
  Pi (0)  e  Pi (  )  e    Po (0)  e  Po (  )  e   0

f bnd  f bnd

1
4
1
1
1

2
2
2
2 
(

E


E
)

(

E


E
)

( o   i ) Eo  Ei
o o
o o 
 i i 2 i i
2
 8
Davis and McCammon, Journal of Computational Chemistry. 1990. 11. 401.
Xiang et al, Journal of Chemical Physics. 2009. submitted.
46
Dielectric boundary force:
Newton’s third law
0.009
Forces
0.006
Third law X
Analytical X
Third law Z
Analytical Z
0.003
0.000
-0.003
-0.006
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
Theta
Xiang et al, Journal of Chemical Physics. 2009. submitted.
47
Acknowledgements
Profs. David Case, Michael Gilson, Hong-Kai Zhao and Zhilin Li
Drs. Jun Wang, Siang Yip
Chuck Tan, Yuhong Tan, Qiang Lu
Qin Cai, MJ Hsieh
Gabe Ozorowski, Seema D’Souza
Morris Chen, Emmanuel Chanco
NIH/GMS
48