.1) VHF 011.
1815,43“ .

 

.n .

.1 ..

Em... ;

6 “$5.“... a.” ...,
s xugflg. .
u . . .

3.9..

S
3’ 5.1“:

4,0. t

.1 5v. . , . q
K Wu .

. . é

. gmfifil.

,. z a. I...

it :1. 313nm”...

.51 n1 .. I‘ .1

A....:1......J......T..., a. 15%

l..‘: ‘3 \‘I I J“
.13! .2. Ilrl. 1
$3....» .

b.
:1 1.. was.“ flu

.
2 , .
.1
.s. a
in...
iifill.

1:5,, -3“: .i 2.
14:21.11»... ignfli
nth... 15.32 g

Id
Id}

..

U
a: 21.11.".
37.1.1
iv)

"an

‘51.

 

)ibtlrw’lo...‘
. Luv... 11.1 in
11390:.“ .
.1...
.1. “rt? . ..
. .. ... . e . 64:314. .. 3.5

a.“ .

u .1 .1? {1 I}...
I. 1115:1131?
it 415:: 1.1.11.1
1 . 1.. . . 1‘ {31711)

1.11

.4 .... : 1 E

 

““338 LIBRARY
_ ., Michigan State
4"“ 2 University

 

 

 

This is to certify that the
dissertation entitled

INTERFACE METHOD AND GREEN'S FUNCTION
BASED POISSON BOLTZMANN EQUATION SOLVER
AND INTERFACE TECHNIQUE BASED MOLECULAR

DYNAMICS

presented by

WEIHUA GENG

has been accepted towards fulfillment
of the requirements for the

Ph. D. degree in Applied Mathematics

 

 

\

/'/4./’":”’/"/‘/

~..’

 

Major Professor’s Signature

Argg [7 2 ,0 0%

Date

MSU is an affirmative-action, equal-opportunity employer

 

-p—-.—.-.-i—.-.—.—-.-.—.—.—.— —.———.

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 KIProj/AccliPres/ClRC/DateDueindd

 

INTERFACE METHOD AND GREEN’S FUNCTION
BASED POISSON BOLTZMANN EQUATION SOLVER
AND INTERFACE TECHNIQUE BASED MOLECULAR

DYNAMICS

By
WEIHUA GENG

A DISSERTATION
Submitted to

Michigan State University
in partial fulfillment of the requirements

for the degree of

DOCTOR OF PHILOSOPHY

Applied Mathematics

2008

ABSTRACT
INTERFACE METHOD AND GREEN ’3 FUNCTION BASED POISSON
BOLTZMANN EQUATION SOLVER AND INTERFACE TECHNIQUE BASED
MOLECULAR DYNAMICS
By

Weihua Geng

This dissertation describes the development of a higher order interface and Green’s
function based method for solving the Poisson Boltzmann (PB) equation with both
geometric singularities and charge singularities. Meanwhile, in the framework of the
implicit solvent model, an interface technique based molecular dynamics method is
developed for biomolecules. Some material presented in this thesis is adopted from
related reprints and preprints.

The author starts in Chapter 1 from a bird’s eye view of the current interface
methods, focusing on matched interface and boundary (MIB) method developed in
Prof. Wei’s group in the past five years.

Chapter 2 reviews the history, status, challenge and application of the PB model.
The numerical difficulties of PB equation solvers in handling discontinuous coefficients
and solutions, geometric and charge singularities are described there. Meanwhile, as
the basis of later chapters, the formulation Of PB problem is provided.

Chapter 3 describes the development of the third generation Of MIB based PB
equation solvers, the MIBPB-III, whose development is based on a Green’s function
formalism and is a natural continuation of the MIBPB-I and MIBPB-II, the previous
two generation of MIBPB with the enforcement of interface flux continuity conditions
to the PB equation and the treatment of geometric singularities respectively. In

the present Green’s function formalism, the MIBPB-III, the charge singularities are

transformed into interface flux jump conditions, which are treated on an equal footing
as geometric singularities in our MIB framework. The MIBPB-III is able tO provide
highly accurate electrostatic potentials at a mesh as coarse as 1.2A for proteins.
Chapter 4 describes the application of the MIB method to the development of the
first PB based molecular dynamics (MD) method that directly admits sharp molecular
surfaces. The classical formulation of electrostatic forces is modified to allow the use
of sharp molecular surfaces. Accurate reaction field forces are obtained by directly
differentiating the electrostatic potential. Dielectric boundary forces are evaluated
at the solvent-solute interface using an accurate Cartesian-grid surface integration
method. The electrostatic forces located at reentrant surfaces are appropriately as-
signed to relevant atoms. Extensive numerical tests are carried out to validate the
accuracy and stability of the present electrostatic force calculation. The resulting
electrostatic forces are compared with those in the literature. The present MIB and
PB based molecular dynamics simulations of biomolecules are demonstrated in conju-
gation with the AMBER package. It has been shown that there is a pressing need to
examine the convergence, stability and reliability of previous PB based MD methods.
Chapter 5 is the summary of author’s thesis contributions and a brief descrip-
tion of important future topics in the PB methods. Emphasis is given to potential

developments in mathematical modeling and computation of biomolecular systems.

Acknowledgments

I have been bearing a heart of appreciation and gratitude during the five-year
pursuit of my PhD degree. The completion of this dissertation is greatly attributed
to numerous supports from many people. I firstly would thank Professors Guowei
Wei, Keith Promislow, Zhengfang Zhou, Changyi Wang, and Andrew Christlieb in
the mathematics department for serving on my dissertation committee, along with
Professor Michael Frazier, who was my committee member, and now is as the head of
Department of Mathematics, University of Tennessee. Their criticism and suggestion
of this dissertation will inspire me in many ways in my future work and careers.

Most importantly, I would like to take this opportunity to thank my thesis ad-
viser, Dr. Guowei Wei, who set me on the right path after I have joined his group.
His systematic guidance and constant push not only work as the major sources of the
completion of this dissertation and its corresponding publication, but also cultivates
me from a student to an independent scientific researcher. His support and encour-
agement in every research related aspect sustain my confidence in entering a highly
competitive academic area. At the moment of graduation, as one of his students, I
feel well prepared to the future and greatly recognized by the society.

I will join Department of Mathematics, University of Michigan, for the three-year
postdoctoral position after obtaining my PhD thesis. I thank Professor Krasny for
his support in recruiting me to this position, and I anticipate a great and fruitful
time with him there. I also would like to thank all the professors who wrote support
letter for my job application. They are Professor William Brown, the undergraduate
director of Math Department, who admitted me to Michigan State University seven
years ago, Professor Leslie Kuhn, my first bio-orientated professor, and Professor
Lijian Yan, my girlfriend’s adviser and instructor Of two of my statistical courses.

iv

The support and help from current and former members of Dr. Wei’s group are
indispensable to this dissertation. I owe a great debt of gratitude tO Dr. Yongcheng
Zhou, who has provided in time suggestions and help to me no matter when he was
a PhD student at MSU or at his postdoctoral position at UC San Diego. He is like a
trailblazer to me in research and career. I am also grateful to the support and help
from Dr. Sining Yu. As the collaborator of several research projects and coauthor of
several papers, I learned many things from her, from the spirit of pursuit of dreams
to the attitude for challenging difficulties. There are many words could be written
for my thanking to all group members for their help and support. To save the place,
I will simply list their names as Professor Shan Zhao, Dr. Yuhui Sun, Dr. Yi Mao,
Mr. Zhan Chen, Mr. Duan Chen and Miss Qiong Zheng.

My PhD degree will not be completed without the help and support from my
family. My father, Zhaoquan Geng, mother, Hong Su, and bother Guangyu Geng
have been giving their unfiagging support from the other side of the earth. They
trust my every decision and appreciate my every improvement. My fiancee, Qiongxia
Song, as the study and life partner with me for my recent two years, is the one I can
really rely on when I have to face difficulties and challenges. Her fantastic cooking
and dedicated care of my living provided endless energy for my work. I owe them a
great deal.

Michigan State University is such a great university to study, research and live. I
appreciate everything I have experienced and everybody I have met during my seven
years stay here. In particular, I am grateful to the continuous financial supports from
the Department of Mathematics in terms of teaching assistantship and from Professor
Wei’s NSF grants in terms of research assistantship. I would like to add a few names

for the people who brought me such great experience here as Xiaoyan Shi, Xijin Yan,

Weixing Song, Lin Liu, Tingting Yi, Jing Wang, Linna Ye, Dongsheng Wu, Adam
Goyt, Karen Brooks, Jaylen Jones, Lei Wang, Likun He, and Xiaoyue Luo with my

best wishes.

vi

Table of Contents

Chapter 1 Introduction to Elliptic Interface Methods 1
1.1 Status and challenges .............................. 1

1.2 Matched interface and boundary method for elliptic equations with discon-
tinuous coefficients and singular sources .................... 4
1.2.1 Introduction ............................... 4
1.2.2 MIB Schemes .............................. 6
Chapter 2 Introduction to the Poisson Boltzmann Equation 13
2.1 Status and challenges .............................. 13
2.2 Formulation ................................... 17

Chapter 3 Treatment of Singularities in the Poisson-Boltzmann Equation 21

3.1 Introduction ................................... 21

3.2 Theory and algorithm .............................. 24

3.2.1 Treatment Of geometric singularities ................. 24

3.2.2 Treatment Of charge singularities ................... 28

3.2.3 The computation of electrostatic free energy of solvation ...... 35

3.3 Results and discussions ............................. 36

3.3.1 Units of the PB equation ........................ 37

3.3.2 Validation ................................ 38

3.3.3 Applications ............................... 46

3.4 Conclusion .................................... 54

Chapter 4 MIB Based Molecular Dynamics Method 65

4.1 Introduction ................................... 65

4.2 Theory and algorithm .............................. 68

4.2.1 Free energy functional ......................... 68

4.2.2 Electrostatic forces ........................... 70

4.2.3 Force expressions involving sharp interfaces ............. 73

4.2.4 The MIB solution of the Poisson Boltzmann equation ........ 75

4.2.5 Dielectric boundary force calculation ................. 81

4.2.6 Special issues about MIBPB-III .................... 90

4.3 Validation .................................... 92

4.3.1 Validation of the surface integration in Cartesian grids ....... 92

4.3.2 Validation of force evaluation ..................... 94

4.4 Application to molecular dynamics simulation ................ 103
4.4.1 The implementation of MIBPB based PB forces in the AMBER

package for molecular dynamics .................... 103

4.4.2 Molecular dynamics simulations .................... 107

4.5 Concluding remarks ............................... 110

Chapter 5 Summary of Thesis Contribution and Future Work 114

vii

5.1 Thesis contribution ...............................

5.2 Future work

References

viii

List Of Tables

3.1

3.2

3.3

3.4

3.5

3.6

3.7

3.8
3.9
3.10

4.1
4.2

Comparison of electrostatic free energies and surface potential errors
for a unit charge in a dielectric sphere (Case 1) ............ 56

Numerical errors of surface electrostatic potentials for an off-centered
charge in a dielectric sphere (Case 2). ................. 57

Component analysis of surface electrostatic potentials for an off-centered

charge in a dielectric sphere (Case 2). ................. 58
Numerical errors of surface electrostatic potentials for two off-centered
charges in a dielectric sphere (Case 3). ................. 59
Numerical errors Of surface electrostatic potentials for six charges in a
dielectric sphere (Case 4). ........................ 59
Errors and convergence of MIBPB-II and MIBPB-III for solving Eq.
(3.11) with molecular surfaces ...................... 60
The E1 errors and convergences Of MIBPB-II and MIBPB—III for solv-

ing Eq. (3.11) with the molecular surfaces of twenty four proteins . . 61
Electrostatic free energies of solvation of twenty four proteins . . . . 62
Electrostatic free energies of solvation Of twenty four proteins ..... 63
Electrostatic potential values and differences (kcal/mol/ec) for cytochrome
C551. ................................... 64
Integration errors for spherical surfaces ................. 112

Comparison of surface areas calculated from numerical integrations and
the MSMS. ................................ 113

ix

List Of Figures

1.1 An illustration of the MIB method. (a) Fictitious domains constructed
with irregular grid points in a 2nd-order MIB scheme; (b) MIB scheme. 8

2.1 Domains of the PB equation (a) Three-domain Model; (b) Two-domain
Model ................................... 19

3.1 Illustration of schemes applied in MIBPB-II (a) an example of geo-
metric singularity; (b)auxiliary points setup in MIBPB-I; (c)auxiliary
points setup in MIBPB-II ........................ 25

3.2 450 extension schemes (a) and illustration (b) .............. 32
3.3 Schemes to find partial derivatives: (a) regular cases; (b) an special case 33

3.4 Comparison of absolute maximum errors (E1) of the surface potentials
of 24 proteins, which are ordered with ascending radii of gyration: 1ajj,
2pde, lvii, 2erl, lcbn, lbor, 1bbl, 1fca, luxc, lshl, 1mbg, lptq, 1vjw,
lfxd, 1r69, 1hpt, lbpi, 451C, 1a2s, 1frd, lsvr, lneq, 1a63 and 1a7m.
Errors of MIBPB-III (green), and PBEQ (blue) at 24 proteins are
respectively linked with dash-lines and solid—lines. Errors of MIBPB—II
(red) are unlinked. Here, D, 0, and A, are respectively for h = 1.0, 0.5
and 0.2521 .................................. 47

3.5 Comparison of electrostatic free energies of solvation of 24 proteins
listed in the order Of gyration radii. (a) Solvation energies; (b) Dif-
ferences of solvation energies between mesh sizes of h = 0.25 and
h. = 0.5A; (c) Differences of solvation energies between MIBPB-III
at h. = 0.25A and the results in the legend; (d) CPU time. ...... 49

3.6 Analysis of surface electrostatic potentials of cytochrome C551. (a)
¢MIBPB-III(h=0.25)§ (b) ¢MIBPB-111(h=o.5) ' ¢MIBPB-III(h=O.25)§ (C) ¢MIBPB-Il(h=0.5)
‘ ¢MIBPB-11(h=0.25); (d) ¢PBEQ(h=0.5) ' ¢PBEQ(h=0.25) --------- 52

3.7 Differences of surface electrostatic potentials of cytochrome C551. (a)
¢MIBPB-III(h=1.0) - ¢MIBPB-111(h=o.25); (b) ¢MIBPB-Il(h=0.5) ' ¢MIBPB-III(h=O.25)§
(C) ¢PBEQ(h=0.25) ' ¢MIBPB-III(h=0.25)i (d) ¢PBEQ(h—_—0.5) ' ¢MIBPB-III(h=o.25)- 53

4.1 An illustration of MIB schemes ..................... 76

4.2 Dielectric boundary force calculation ................... 87

4.3

4.4

4.5
4.6
4.7
4.8

4.9

4.10

4.11

Comparison of electrostatic forces of a two—atom systems Obtained at
h = 0.21A. Red solid line: present MIB based method with molecular
surfaces; Blue dotted line: PBEQ with smooth interface. (a) Solvation
free energy of the system; (b) Total force on the moving atom calcu-
lated by using solvation free energies; (c) Total force on the moving
atom calculated from the summation of reaction force and dielectric
boundary force; ((1) Reaction force; (e) Dielectric boundary forces; (f)
Dielectric boundary forces when the moving atom on x-axis is in be-
tween 3.021 and 4.021 ...........................

Electrostatic forces calculated for two separating charged atoms. The
vertical axis is the position of the second atom and the horizontal is the
forces. (a) PBEQ at h = 0.2121; (c) PBEQ at h = 0.4211; (b) MIBPB
at h = 0.21A; (d) MIBPB at h = 0.4221. ................

Forces calculated for a three-atom system ................
An eighteen—atom molecule Cyclohexane ................
Homo sapiens ...............................

Correlations of forces obtained at different mesh sizes by the MIBPB
for protein 1aj j. (a) Reaction field forces between mesh sizes of 0.25 and
0.5A; (b) Reaction field forces between mesh sizes of 0.25 and 1.0A; (c)
Toric forces between mesh sizes of 0.25 and 0.5A; (d) reentrant forces
between mesh sizes of 0.25 and 0.5A; (e) Dielectric boundary forces
between mesh sizes of 0.25 and 0.5A. ..................

Correlations of forces Obtained at different mesh sizes by the PBEQ for
protein 1ajj. (a) Reaction field forces between mesh sizes of 0.25 and
0.5A; (b) Reaction field forces between mesh sizes of 0.25 and 1.0A;
(0) Dielectric boundary forces between mesh sizes Of 0.25 and 0.5A. .

A comparison of 10 ps MD simulations of the alaning dipeptide ob-
tained with PB based methods. Results of the MIBPB and AMBER
PB are illustrated in red and blue lines, respectively. (a) RMSD; (b)
Solvation free energies; (c) Total potential energies; (d) Temperature.

A comparison of 10 ps MD simulations of the WW domain binding
protein-1 obtained with PB based methods. Results of the MIBPB
and AMBER PB are illustrated in red and blue lines, respectively.
(a) RMSD; (b) Solvation free energy; (c) Total potential energy; (d)
Temperature ................................

xi

96

99
100
102
102

104

105

108

109

Chapter 1

Introduction to Elliptic Interface

Methods

1.1 Status and challenges

The numerical solution of elliptic equations with discontinuous coefficients and sin-
gular sources has drawn much attention in recent years [12, 13, 20, 29, 35, 44, 62,
63, 66, 94, 110, 116]. Simple Cartesian grids are preferred in these studies since the
complicated procedure of generating unstructured grid could be bypassed, and well
developed fast algebraic solvers could be utilized. The importance of elliptic interface
problems has been well recognized in a variety of disciplines, such as fluid dynamics
[31, 39, 60, 91], electromagnetics [50, 51] and material science [57]. However, to con-
struct highly efficient methods for these problems is a difficult task due to the low
global regularity of the solution. Traditional numerical methods that are constructed
with the assumption of smooth solutions cannot perform at designed accuracy, and
might even diverge.

A pioneer work in this field was due to Peskin, who proposed the immersed bound-

ary method (IBM) in the 1970—19808 [94—96] to model complicated time-varying ge-
ometric boundaries via a singular source on the interfaces. The success of the IBM
method is due to its flexibility, efficiency and robustness. Tornberg and Engquist [110]
extended Peskin’s method by using a high order approximation to the delta function.
However, the IBM approaches typically smear the interface.

Many high order boundary schemes have been developed. Gibou and Fedkiw [44]
constructed a fourth-order boundary scheme to solve the Laplace and heat equations.
Linnick [80] and Fasel developed another fourth-order boundary scheme for flow past
cylinders. These high order schemes have a potential to be generalized for elliptic
problems with discontinuous coefficients.

The immersed interface method (IIM), proposed by LeVeque and Li [73] is a
remarkable second order sharp interface scheme. The IIM has been regarded as a
major advance in the field. It has been made robust and efficient over the past
decade [1, 28, 74, 75, 104]. The IIM formulations in polar coordinates [78] and using
a finite element method [76] were developed. A fast IIM algorithm was constructed
[74] for interface problems with piecewise constant coefficients.

The ghost fluid method (GF M) [32], proposed by Osher and his coworkers, is a
relatively simple and easy-to—use approach. It is developed to treat contact disconti-
nuities in the inviscid Euler equations. The GFM is typically first-order accurate for
interface problems, including the elliptic one [82] and could be of second-order accu-
racy for elliptic irregular domain problems. It directly incorporates jump conditions
into the numerical discretization such that the symmetry of the FD coefficient matrix
is maintained, allowing the use of standard fast solvers. In high dimension, the jump
in the normal derivative is correctly captured but jumps in the tangential derivatives
are neglected [82] so that the GFM can be applied dimension by dimension.

2

The finite element method, proposed by Babuéka [3] and many other researchers
[17, 76], is a nature approach to construct a solution for irregular interfaces. In
particular, a discontinuous Galerkin based technique is proposed by Guyomarch [49].
A relevant, while quite distinct approach is the integral equation method for complex
geometry, proposed by Mayo, Greengard and coworkers [88, 89]. In this approach, the
ghost cell is introduced outside the domain as a fictitious domain. Aforementioned
methods have found much success in scientific and engineering applications [2, 12, 13,
29, 35, 44, 56—59, 62, 66, 72, 77, 78, 80, 83, 103, 104, 113, 116, 119].

One of the most challenging problems in the field is the solution of elliptic equa-
tions with sharp-edged interface, i.e., non-smooth interfaces. Numerical solutions to
this class of problems have widespread applications in science and engineering, such
as electromagnetic wave scattering and propagation [21, 40, 93], wave-guides analysis
[90], plasma-surface interaction [87], friction modeling [111], and turbulent-flow [5].
To the best of our knowledge, none of the aforementioned methods proposed for e1-
liptic interface problems have been directly applied to the treatment of sharp-edged
interfaces. Essentially, as the gradient near the tips Of sharpedged interface is not
well defined, some earlier interface methods might not work. Most existing results on
this class of problems are obtained by using finite element methods [90, 93]. However,
finite element methods might exhibit a reduced convergence rate when used for the
analysis of geometries containing sharp edges [56, 93]. Consequently, dramatic local
mesh refinement is required in the vicinity of sharp edges [22], and leads to severe
increase in computational time and memory requirement. In particular, local mesh
refinement does not work if the solution is highly oscillatory due to the so called pol-
lution effect [4], which is a common situation in dealing with electromagnetic wave
scattering and propagation. Hou and Liu proposed a finite element formulation [56]

3

for solving elliptic equations with sharp-edged interfaces. Remarkably, these authors
have achieved about 0.8th order convergence with non-body-fitting grids.

Another challenging problem is to realize the higher-order accuracy in high dimen-
sion, which is particularly desirable for problems involving both material interfaces
and high frequency waves where conventional local adaptive refinement approaches
do not work well. Typical examples are the interaction of turbulence and shock, and
high frequency wave propagation in inhomogeneous media [10]. Simple Cartesian
grids are preferred in these situations because of the bypass of the mesh generation,

better temporal stability and the availability Of fast algebraic solvers.

1.2 Matched interface and boundary method for
elliptic equations with discontinuous coefficients

and singular sources

1 .2.1 Introduction

The development of MIB schemes in solving elliptic problems with singular inter-
faces are described as the follows. It was first proposed by Zhao and Wei [126] as a
systematic higher-order method for electromagnetic wave propagation and scattering
in dielectric media. After that, for straight interfaces, MIB schemes of up to 16th
order have been constructed [126, 129]. For lightly curved interfaces, up to 6th or-
der schemes have been demonstrated [129]. Meanwhile, it has been generalized for
solving elliptic equations with curved smooth interfaces by Zhou et al [130]. The

MIB approach makes use of fictitious domains so that the standard high order cen-

tral finite difference (FD) method can be applied across the interface without the
loss of accuracy. However, the earlier MIB scheme does not maintain the designed
accuracy due to the presence of geometric singularities or large curvature [128-130].
It is very difficult to obtain the fictitious values around geometric singularities and
interfaces with large curvature because the number of unknowns around these areas
is usually greater than the number of interface jump conditions. Recently, the afore-
mentioned challenge is overcome by making use of additional set of jump conditions
and iteratively solve fictitious values around geometric singularities [122, 123]. It was
demonstrated that 3D 4th—order MIB schemes can be constructed for interfaces with
moderate geometric singularities. Such schemes are particularly efficient for prob—
lems involving both material interfaces and high frequency waves. To speed up the
MIB schemes, a systematic strategy is developed to Optimize the structure of the
MIB matrix, which dramatically accelerates the convergence rate of the linear system
[122]. The representations of fictitious values are improved by selecting different set
of points around irregular points. This improvement makes the MIB matrix opti-
mally symmetric and banded, and consequently reduces the computational time. It
is shown that the optimized MIB method converges fast even for large systems where
the earlier MIB method converges 100 times slower or even diverges.

Parallel to the development of its schemes, the MIB method has been applied to
calculate the electrostatic potential and free energy of biomolecules. The interface
problem in the Poisson-Boltzmann equation was first addressed by an MIB based
(MIBPB-I) solver to explicitly consider the flux continuity condition in the finite dif-
ference framework [128]. However, the earlier MIBPB-I solver does not maintain the
2nd order convergence for molecular surfaces of proteins because of the existence of
cusps, sharp edges, sharp wedges or self-intersecting surfaces. Moreover, due to the

5

asymmetric matrix of the interface method, the MIBPB-I requires a large number of
iterations in solving linear equations. In fact, its matrix does not converge for large
proteins or small proteins with dense grids. Another MIB approach, called MIBPB-II
[121], overcomes these difficulties. Apart from its ability to maintain the designed
2nd order convergence under the presence of geometric singularities, the MIBPB-II
has an optimally symmetrical matrix, which dramatically reduces the number of it-
erations. However, when the mesh size approaches half of the van der Waals radius,
the MIBPB-II cannot maintain its accuracy because the grid points that carry the
interface information overlap with those that carry distributed singular charges. A
Green’s function method is developed to tackle this problem in this thesis. In the
present Green’s function formalism, the charge singularities are transformed into in-
terface flux jump conditions, which are treated on an equal footing as the geometric
singularities in our MIB framework. The resulting method, denoted as MIBPB-III, is
able to provide highly accurate electrostatic potentials at a mesh as coarse as 1.2A for
proteins (i.e., the smallest atomic radius in protein). Consequently, at a given level of
accuracy, the MIBPB-III is about three times faster than the APBS, a recent multi-
grid PB solver. Detailed argument related to MIBPB solvers will be given in Chapter

3.

1.2.2 MIB Schemes

To illustrate the schemes of MIB method, let us consider an elliptic equation of the

form

-V'(€(F)V¢(rl) + N2(r)¢(r)=f(r)- (1-1)

The interface F divides the whole domain into two separated parts, Q“ and (2+. The

jump conditions across the interface are defined as

[4] WM - 45—(1‘), (1-2)

6+(r)V¢+(r) - n — e- (r)V¢— (r) - n (1.3)

l€¢nl

where n = (n$,ny,nz) is the outer normal direction of the interface I‘. Here [(25]
and [64511] are given or can be computed from other quantities. Appropriate far-field
boundary conditions are also prescribed and the outer boundary of 9+ is chosen to
be regular. A standard Cartesian grid is used along with the standard finite differ-
ence schemes. However, near the interface, the standard finite difference schemes lose
their designed convergence. We therefore implement the interface jump conditions to
restore the accuracy. To this end, we classify all the grid points into two classes, the
regular ones and the irregular ones. An irregular grid point is defined as one where
the standard finite difference scheme involves grid points across the interface. Obvi-
ously, for a given grid point near the interface, it may and may not be an irregular
point, depending on the order of the finite difference scheme used. For a given finite
difference scheme, all the irregular grid points on one side of the interface constitute
a fictitious domain, see Fig. 1.1(a). Fictitious values, which are the linear combi-
nation of regular points and jump conditions, at irregular points are computed as a
continuation of the solution from the other side of the interface by using the jump
conditions. The fictitious value is used when a finite difference scheme reaches across
the interface. In the follows, we discuss how to calculate fictitious values.

Assume that the interface I‘ intersects the grid in the :r-direction at a point

(2'0, j, k), which is located between (2', j, k) and (z' + 1, j, k). We therefore have two

   

i i+1 i+2 1+3 I+4

a (b)
Figure 1.1: An illustration of the MIB method. (a) Fictitious domains constructed
with irregular grid points in a 2nd-order MIB scheme; (b) MIB scheme.

irregular grid points, (2', j, k) and (i + 1, j, k), for the second order finite difierence
scheme in the :r-direction near the interface. The fictitious values at these two irreg-
ular points, f+(i,j, k) and f‘ (i + 1,j, k), are to be determined. However, one of the
jump conditions, Eq. (1.3), is defined in the normal direction of the interface at point
(i0, j, k). It is convenient to introduce a local coordinate (5, 11,0 such that { is along
the normal direction, and n is in the a: - y plane. The coordinate transformation can

be given as

f 1:
77 = P‘ y 1 (1-4)
( z

where p is the transformation matrix

sin w cos 6 sin 11) sin 0 cos 't/}
P = —sin0 0089 0 - (1-5)

—cos¢cos€ —cost/)sin6 sintb

I- d

 

 

Here 6 and 1,!) are the azimuth and zenith angles with respect to the normal direction
5, respectively. By differentiating Eq. (1.2) along two tangential directions, 7) and C,

we can generate two additional jump conditions

[051,] = ($3121 + $31022 + 625:1923) - (49.351921 + 05,71022 + 452—1223) (1-5)

W] = (@181 + $31032 + ¢>L~+P33l - (@331 + $51732 + 4.2—P33). (1.7)

where p0. is the z'jth component of the transformation matrix p. Since Eq. (1.4)

implies

 

 

 

 

95g 493:
$77 = p ' 96y 7 (18)
9% C52

Eqs. (1.3), (1.6) and (1.7) can be rewritten as follows:

 

 

 

 

r):
‘ ‘ (PE
lcétl ,+
19y
[4572] =C' ’ (19)
ab;
95
[cl <25?
P;
where
7 T '- 7
C1 P116+ —P11€_ P126+ -P12t— P13€+ -P13€_
C2 C2 2 P21 -P21 P22 ~P22 P23 -P23 ’ (1'10)
C3 L P31 -P31 P32 -P32 P33 —P33

 

 

 

 

where C,- represents the 2th row of matrix C. It is clear that none of these three jump
conditions is easy to implement because there are six derivatives, (15;, qt; , 92?, qt; , (1):,
and a; to be calculated in appropriate subdomains near the interface. For complex
solvent-molecule interfaces of macromolecules, it is often very difficult to evaluate
some of these derivatives. Therefore, in the MIB method, we avoid calculating two
most difficult derivatives by eliminating them with two relevant jump conditions. Fig.
1.1(b) illustrates a case where the <ng is difficult to compute and is to be eliminated.

In general, after the elimination Of the lth and mth elements of the array (Oi, 9'5; ,

10

05+, $37, 05:, 452—), Eq. (1.9) becomes

a [gag] +b [an] +c [44] = (aC1 +bC2 +cC3)- , (1.11)

 

 

where

a = C2IC3m - C3IC2m (H?)
b = C3ZCIm “ 01103171

C : CIIC2m " CZICIm-

We therefore use Eqs. (1.2) and (1.11) to determine two fictitious values near the in-
terface along a specific mesh line at a time. This procedure is systematically repeated
to determine fictitious values along other mesh lines and at other interface locations.
In this manner, we have effectively reduced a 3D interface problem into a 1D-like one.

Returning to the situation in Fig. 1.1(b), after eliminating ab; , (f); is evaluated
by using interpolation schemes with solution values from the neighboring points. In
details, (153’ at (10,], k) will be interpolated by using grid value at points (imj —- 1, k),
(2'0, j, k) and (to, j + 1, k.) Since these three points are not normal grid points, we need
to further interpolate the values at (to, j — 1, k) by using grid values at (z + 1, j — 1, k),

(1+ 2,j —- 1, k) and (1+ 3,j — 1, k), and the value at (10,)" + 1, k) by using grid values

11

at (2+ 1,j + 1, k), (2' + 2,j + 1, k) and (2 + 3,j + 1, k). For the value at (2'0,j, k), the
fictitious value at (2', j, k) and grid values at (2+ 1,j, k) and (2 +2, j, k) will be applied
for the interpolation. For a given local geometry, we will choose to compute one of
two partial derivatives, (1'): and (15;, such that the evaluation in the a: — 2 plane can be
easily carried out. In fact, care is required to select appropriate interpolation schemes
so that the resulting MIB matrix is optimally symmetric and diagonally dominant
for arbitrarily complex solvent-r'nolecule interfaces with geometric singularities. A
detailed discussion of the matrix optimization and 3D MIB methods for geometric
singularities is given in Ref. [122]. A rigorous validation of MIB based PB solver,

MIBPB-II, for molecular surface singularities can be found in Ref. [121].

12

Chapter 2

Introduction to the Poisson

Boltzmann Equation

2.1 Status and challenges

Electrostatic interactions are omnipresent in biomolecules [69, 112]. About 45% of
aminoacids in globulare proteins are either ionized under physiological conditions or
with polar groups in their sidechains [36]. A high electron charge density Of one
charge per 1.7A is found in double-strained DNA chains. Electrostatics plays a
paramount role in the structure, function, stability, and dynamics of macromolecules,
such as signal transduction and other metabolic processes. The interaction of each
pair of (partially) charged particles is governed by the Coulomb potential, which
has a non-vanishing impact over a wide range of distance. Therefore, accurate and
efficient evaluation of electrostatics has been one of the most challenging issues in
computational molecular biology. Although a variety of methods, including Ewald
summations, Euler summations, periodic images and reaction field theory, have been

developed in the past few decades, explicit electrostatic calculations of biomolecules

13

in the solvent remain extremely expensive. Alternative implicit solvent models [6, 33]
have become very popular [26, 38, 55, 101, 105, 114, 124] since the pioneering work
by VVarwicker and Watson in the early 19805 [118], and Honig in the 19908 [55]. The
implicit solvent theory retains a microsc0pic treatment of biomolecules, while adopts a
macroscopic mean-field description of the solvent. Compared with the explicit model,
implicit model sacrifices some molecular details Of the solvent but greatly reduced the
computational demand. Moreover, for some macroscopic average solvent properties
such as dielectric effects, that are difficult to simulate accurately with explicit solvent,
can be built into an implicit solvent model [97].

In such an approach, the Poisson-Boltzmann equation, or the Poisson equation
if no salt is present, is solved for electrostatic potentials. It provides a continuum
model at the equilibrium state to describe solvent by dielectric constant and the mo-
bile ions by the Boltzmann distribution Of its ionic density. The PB equation is the
mathematical basis for the Gouy-Chapman double layer (interfacial) theory; first pro—
posed by Gouy in 1910 [47] and complemented by Chapman in 1913 [23]. It has also
proved to be a successful quantitative tool yielding accurate descriptions of electri-
cal potentials, diffusion limited processes, pH-dependent properties of proteins, ionic
strength-dependent phenomena [55]. The solvation free energies of organic molecules
[37, 70, 107], pKa values [43, 92, 117], and electrostatic forces for molecular dynamics
[85, 97] can be calculated via the solution to the PB equation.

Mathematically, the PB equation is an elliptic equation with discontinues coef-
ficients and singular source terms that describes electrostatic interactions between
molecules in ionic solutions. Although the PB equation can be analytically solved
for a few simple cases, it relies on numerical approaches to obtain useful solutions
for realistic biological systems. A vast variety of computational approaches, such as

14

finite difference methods (F DM) [38, 55, 61, 84—86, 105, 118], finite element methods
(F EM)[8, 9, 26, 53], and boundary integral methods (BIM) [14—16, 19, 65, 79, 114,
124], have been developed in the past few decades.

All methods are subject to certain inherent advantages and limitations, which
are closely tied to the associated underlying formulation. FEMS are optimal for
applications that require the rapid adaptation of grid points to account for struc-
tural variation of biomolecules. However, generating unstructured grids for complex
biomolecular interfaces is time-consuming, especially for large biomolecules. BIMs
enjoy several intrinsic advantages such as fewer unknowns, exact far-field treatment,
and accurate representation Of surface geometry and charge singularity. In light Of
these advantages, the BIM, especially when accelerated with the fast methods, can
also provide an efficient computational approach. However, BIM is not very efficient
in dealing with the nonlinear term in the PB model. FDMS are the main workhorse
for solving the PB equation in computational structural biology and computational
genetic engineering for the following reasons. ( 1) Using 3D Cartesian grids, FDMS
avoid the time-consuming grid generation step. (2) 3D Cartesian grids are a stan-
dard option in most programming codes; therefore, finite-difference based PB solvers
naturally fit into these simulation packages. (3) Compared to the BIM, the solution
is obtained directly on the grid points of interest, and is readily available for being
used in computing the electrostatic forces or solvation energy, (4) FDMS are rela-
tively simple, and particularly in conjunction with multigrid linear algebraic solvers,
can offer the best combination of speed, accuracy, and efficiency, making them the
most popular approach [8]. In the past decade, the main developments with respect
to the PB methodology have been focused on acceleration of these numerical methods
for the PB equation. It was argued by Baker that finite-difference based PB solvers,

15

particularly in conjunction with multigrid linear algebraic solvers, can offer the best
combination of speed, accuracy, and efficiency, making them the most popular ap-
proaches in structural biology [6]. Among these approaches, Cartesian grid based
finite difference methods are commonly used in many popular software packages,
such as DelPhi [68, 100], UHBD [27], MEAD [100], APBS [7, 52, 53] and CHARMM
[18].

Implicit solvent models require a solvent-molecule interface to separate the sol-
vent domain from the biomolecular domain. The molecular surface (MS) [25, 71]
is commonly used in the PB equation for this purpose. The dielectric constants
of the solvent domain and molecular domain are usually chosen as 80 and 1 (or
2), respectively, leading to discontinuous coefficients in the PB equation. Moreover,
molecular surfaces admit geometric singularities [25, 30, 46, 102], such as cusps and
self-intersecting surfaces. Explicit interface treatment of geometric singularities has
not been considered in the literature. Consequently, no Poisson-Boltzmann (PB)
solver of 2nd-order convergence, i.e., the numerical error reduces by a factor Of four
when the grid spacing is halved, has ever been reported in the context of molecular
surfaces of macromolecules. Indeed, irregular interfaces and geometric singularities
cause numerical instability and slow down the convergence of existing PB solvers.
Therefore, enhancing the stability and accelerating the convergence of PB solvers are
pressing mathematical issues in developing the next generation PB solvers.

In summary, as a well established mathematical mode with wide application in
biophysics, biochemistry and biology, The PB equation encounters numerical and
computational challenges. The research conducted in this thesis focuses on the con-
struction of Green’s formulation to regularize the charge singularities, which co—exist
with geometric singularities, and the development of MIB based molecular dynamics

16

methods.

2.2 Formulation

Consider an open domain I? E R3. Let P12 and F23 be two disjoint interfaces which
divide S) into three disjoint open subdomains, Q = Q" UQOUQ+, see Fig. 2.1(a).
Here {2‘ is the biomolecular domain, {20 is a thin ion-exclusion layer, and {2+ is the
solvent domain. By assuming the Boltzmann distribution for the equilibrium ionic
electrostatic potential energy eccb(r) in the dielectric continuum treatment for the

ions in the solvent domain, the Poisson-Boltzmann equation, takes the form [54]

 

ec

. Nm
_V. (c(r)V¢)(r)) + 52(1) (L81) sinh (82:21?) = 47r§q,6(r — r,), (2.1)

where ¢(r) is the electrostatic potential and q, is the (fractional) charge at position rj.
Here, constants ec, k B, T and Nm are respectively the electronic charge, Boltzmann
constant, the absolute temperature and number of charges in the biomolecule. The

dielectric coefficient C(I‘) and the ionic strength k(r), as functions of the space, are

defined as
6‘ r E Q‘;
e(r) : (2.2)
6+ 1‘ E 00 UQ+.
and
0 r E 9" U 520;
K.(r) = (2.3)
R r 6 (2+.

Eq. (2.1) is subject to the following far-field boundary condition

C(00) = 0. (2.4)

17

However, in practical computations, the Direchlet boundary condition

N
3(1) - V: ——qi—e—"‘~lr-ril (2 5)
' —i—1€+lr_ri| ’ '

is used for the linearized PB equation. Due to the discontinuous nature of coefficients
6 and k, Eq. (2.1) should be solved with following additional interface jump conditions

across F12

[Plrlz = 0 (2.6)
[€¢n]r12 = O (2.7)
and across [‘23
l¢lF23 = O (28)
l¢an23 = O (29)

where n is the outer normal direction of the interface. Mathematically, these jump
conditions are required to ensure the uniqueness of the solution. Physically, these
conditions are required by the continuity of the electrostatic potential and its flux
across the interface.

Since the ion-exclusion layer is typically about 1 to 2 A in thickness, the thin
layer domain 90 normally has very limited impact to the electrostatic potential. For
simplicity, we adopt a two-domain model as illustrated in Fig. 2.1(b), where the entire

domain is divided by the boundary Of the molecule, I‘, into two disjoint domains 9‘

18

 

 

\ \
Solvent + + ...———-—-—'Mobile Ions Solvent + + /Mobile Ions

 

Ion-exclusion layer +

 

 

 

 

 

(a) (b)
Figure 2.1: Domains of the PB equation (a) Three-domain Model; (b) Two—domain
Model

and 9+. Correspondingly, we consider the simplified 6(1‘), and rc(r) as

c“ r E 9';
6+ r 6 {2+
and
0 r E Q“,
k(r) = . (2.11)
k r 6 12+

The interface jump conditions across F are given as

l¢lr = 0 (2-12)

[645an = 0- (2-13)

Numerically, it is convenient to solve the PB equation in its dimensionless form

obtained by defining ¢(r) = gig—(152. The resulting dimensionless PB equation is given

19

Nm
—V - (c(r)V¢(r)) + 222(r) sinh(<p'(r)) = Z 2,6(r — r2), (2.14)
i=1

where z, = (i—Ef) q,- is the charge fraction at position ri. Under a weak ionic field

assumption it is convenient to solve the linearized PB equation

Nm

-V ' (€(I‘)V¢(I‘)) + r62(1‘)¢(1F) = 2 225(1‘ - ril- (2-15)

221

In our work, numerical methods are discussed for solving this linearized PB equation.

20

Chapter 3

Treatment Of Singularities in the

Poisson-Boltzmann Equation

3.1 Introduction

The rigorous treatment of discontinuous coefficients and singular sources in elliptic
equations is a challenging task in applied mathematics and scientific computing. Since
Peskin’s pioneer work on immersed boundary method (IBM) [94], a number of other
elegant methods have been constructed. Among these approaches, the ghost fluid
method (GFM) proposed by Fedkiw, Osher and coworkers [32, 82] is relatively simple
and easy to use for complex geometry. Finite element formulations are more suitable
for weak solutions [3, 17, 76, 81]. An upwinding embedded boundary method was
proposed by Cai and Deng [20] for electromagnetic waves in dielectric media. Gibou
and Fedkiw [44] proposed a 4th order extrapolation scheme for solving Laplace and
heat equations on irregular domains. An integral equation approach was developed
for complex geometries [88]. LeVeque and Li [73] proposed a remarkable second

order sharp interface scheme, the immersed interface method (IIM). Their method

21

has been applied to many practical problems, including the 2D PB equation [98].
Recently, a matched interface and boundary (MIB) method has been developed for
interface problems [126, 129, 130]. The MIB is of arbitrarily high order in principle.
Fourth and sixth order MIB schemes have been demonstrated for curved interfaces
[129, 130]. This technique was used to developed our first generation MIB based PB
solver, denoted as MIBPB—I, for arbitrary biomolecular interfaces [128]. Nevertheless,
none Of the aforementioned interface techniques is able to maintain its the designed
order of convergence and accuracy in the presence of geometric singularities. Indeed,
technically, constructing higher order convergence schemes for elliptic equations with
geometric singularities is extremely challenging, despite Of great desire for doing so
in practical applications. To our knowledge, the best result in the literature is of
0.8th order convergence reported by Hou and Liu [56], achieved with a finite element
formulation in 2D. More recently, the MIB method has been extended for solving
elliptic equations with geometric singularities in 2D [123]. Most recently, we have
constructed a 2nd order MIB scheme for solving the PB equation with geometric
singularities of molecular surfaces [121]. Extensive numerical experiments indicate
that our MIB based PB solver, denoted as MIBPB-II, outperforms traditional PB
solvers in terms of accuracy and efficiency [121].

A remaining issue in our MIB based PB solver is the efficient treatment of charge
singularities. Technically, in most finite difference based PB solvers, Dirac delta
functions in the PB equation, or singular charges in a biomolecule, are redistributed
to their neighboring grid points. This treatment works well in commonly used PB
solvers where the major error is produced due to the neglect of interface continuity
conditions. However, it affects the numerical accuracy Of our MIBPB-II solver because
of the possible interference of geometric interface and charge singularities. Stated

22

differently, a set of auxiliary grid points that carry interface jump conditions may also
carry singular charges, which leads to an accuracy reduction. This typically happens
when the mesh size approaches half of the van der Waals radius. Consequently, at a
given level of accuracy, our MIBPB-II delivers highly accurate electrostatic potentials
when the grid spacing is smaller than 0.6A, but it behaves like a normal PB solver
when a coarser mesh is used.

The objective of the present work is to overcome this obstacle by introducing
Green’s function technique into our MIB formulation. The separation of the PB
solution into regular and singular parts was suggested by Zhou et al [131]. The
feedback of the singular part to the regular part was not included in their formulation.
Chern et al have formulated the problem in a more rigorous manner [24]. However,
they have only considered 2D problems with smooth interfaces. In the present work,
we consider a Green’s function formulation for the molecular surfaces of proteins with
possible geometric singularities. The Green’s function approach effectively transfers
the contribution of charge singularities into a set of interface jump conditions which
are determined by solving the corresponding Laplace equation with given singular
sources. The charge-induced interface jump conditions are perfectly compatible with
our MIB techniques for complex interfaces [121, 122]. Therefore, we are able to treat
geometric and charge singularities on an equal footing. The resulting PB solver,
denoted as MIBPB-III, is capable Of eliminating the interference of the geometric
and charge singularities in solving the PB equation, and is of 2nd order convergence.
The present MIBPB-III is able to deliver high numerical accuracy at a grid spacing
as large as the van der Waals radius (about 1.2A) for macromolecules. Consequently,
it is about 3 times faster than APBS [7, 52, 53], a recent multigrid PB solver.

The rest of this chapter is organized as follows. Section 3.2 is devoted to the

23

theoretical formulation and computational algorithm. The treatment Of geometric
singularities is very briefly discussed. A more detailed description of this issue can be
found in Refs. [121, 122] or Dr Sining Yu’s PhD thesis. The Green’s function formu-
lation of charge singularities is developed in the framework of our MIB method. A
second order MIB scheme is constructed to solve the charge-induced boundary value
problem with molecular surface singularities. The resulting MIBPB-III is extensively
validated in Section 3.3. The benchmark tests, such the Kirkwood sphere [67], the
molecular surfaces of few-body systems, and twenty four proteins are used to examine
the accuracy and test the speed of convergence of the proposed MIBPB-III. Compar-
isons are given to our previous MIBPB-II and two other established PB solvers. This

chapter ends with a brief conclusion summarizing the main points.

3.2 Theory and algorithm

3.2.1 Treatment of geometric singularities

Implicit solvent models require a biomolecular surface definition and a prescription of
dielectric functions both inside and outside the biomolecdle. The solution of the PB
equation is very sensitive to the discontinuous dielectric function. Earlier biomolec-
ular surfaces, such as van der Waals surfaces solvent accessible surfaces [71], are
not smooth and cause instability in numerical simulations. The molecular surface
(MS) as proposed by Lee and Richards [71] are designed to provide smoother inter-
faces. In fact, most PB calculations have used the MS. However, in certain geometric
situations, the MS definition admits cusp and self-intersecting surface singularities

[25, 30, 46, 102], which still cause numerical instability. Smoothly varying functions

24

 

I 1+1

Figure 3.1: Illustréffon of schemes applied in MIBPB-II (a) an example of geo-
metric singularity; (b)auxiliary points setup in MIBPB-I; (c)auxiliary points setup in
MIBPB-II

were proposed to avoid this problem [48, 61, 85]. Bimolecular surfaces are some of the
most important concepts in molecular biology. The stability and solubility of macro-
molecules, such as protein, DNA and RNA, are determined by how their surfaces
interact with solvent and/or other surrounding molecules. Therefore, the structure
and function of macromolecules depend on the features of their biomolecular surfaces.
Nonetheless, McCammon and coworkers [108] have shown that some atomic centered
dielectric functions may lead to unphysical interstitial high dielectric regions in im~
plicit solvent models. The solvent exclusion property of the MS is crucial to implicit

solvent models.

25

The MIB based PB equation solver, MIBPB-I, was the first Of its kind that rig-
orously treats the flux continuity conditions in bimolecular context. The MIBPB-II
was designed to maintain the second order convergence even if the molecular sur-
face admits geometric singularities. Here we will briefly describe the improvement
of MIBPB-II over MIBPB-I in their corresponding schemes. In MIB method, we
need to find the fictitious values f + and f ‘, which are the linear combination Of
neighboring points, as illustrated in Section 1.2.2. However, with the increase of the
sharpness and singularities of the interface, the following difficulties will be encoun-
tered, resulting in the failure of Obtaining fictitious values f + and f ". First, the
point, where its fictitious value needs to be Obtained is isolated by two points in the
other domain as illustrated in Fig. 3.1(a). Second, not enough auxiliary points are
obtainable from the desired domain. The details of solving these difficulties can be
found in [122, 123]. Briefly, the MIBPB-II is designed to handle these difficulties in
the following sequences and treatment. First, the scheme checks if the fictitious value
can be obtained from the other directions. Second, the scheme combines the four
jump conditions on both sides of the point e.g. the jump conditions on 232,313,275, and
.730, ya, 20 in Fig. 3.1(a) and the points on which the fictitious values are obtainable
to find the fictitious values at the points with difficulties. With these two treatments,
most singular geometry can be handled.

Another progress in MIBPB-II is at the setup of the six auxiliary points to speed
up the convergence. Fig. 3.1(b) shows the points setup in MIBPB-I and Fig. 3.1(C)
gives that in MIBPB-II. This simple alternation makes the final linear matrix having
a smaller condition number and better symmetry, which improves the convergence
rate in great scale [121, 122].

The accuracy and convergence order of MIBPB-II in the presence of geometric

26

singularities are carefully tested and validated by a few test cases similar to what was
used in next sections. At the same time, the speed of convergence in terms of the
number of iterations NBiCG is also studied. We will just briefly describe the results as
the major work in this chapter is the treatment of charge singularities. In [121], the
MIBPB-II is applied to electrostatic potential calculations of a dielectric sphere with a
unit charge, which admits an analytical solution [67], and a set of twenty four proteins.
As a comparison, the finite difference based PBEQ and APBS are employed for the
same electrostatic calculations. For the analytical case, the MIBPB—II at a coarse grid
of h = 0.5A was found to be more accurate than the PBEQ and APBS at fine grid
of h = 0.05A in terms Of both surface electrostatic potential and solvation energy.
For twenty four proteins, we found that the solvation free energies computed by the
MIBPB-II has essentially converged at the grid resolution of h = 0.5A. The results
of APBS and PBEQ converge toward those of the MIBPB-II when their grids are
refined. The discrepancies between converged MIBPB—II solvation free energies, which
are obtained at the grid of h = 0.5A, and those of PBEQ obtained at h = 0.15A are
nearly 20 kcal/mol for proteins with a large radius of gyration. Slightly smaller
discrepancies are found between converged MIBPB-II solvation free energies and those
of APBS obtained at about h ~ 0.2A. The accuracy of both PBEQ and APBS shows
a dependence on the radius of gyration. The present MIBPB-II requires more CPU
time than APBS and PBEQ at given grid resolution, but less CPU time at a given
accuracy.

As the interface flux jump conditions are rigorously enforced at complex solvent-
solute interfaces and convergence order is systematically restored at geometric singu-
larities, the proposed MIBPB-II solver ought to deliver more reliable surface electro-
static potentials. Indeed, converged surface electrostatic potentials are obtained at

27

the resolution of h = 0.5A. The discrepancies between converged MIBPB-II surface
electrostatic potentials and those of PBEQ at the same grid resolution are about
5kcal/mol/ec. Such discrepancies would induce a consequence in a quantitative anal-

ysis.

3.2.2 Treatment Of charge singularities
Green’s function formulation of the Poisson—Boltzmann equation

N umerically, second and higher-order numerical implementation of Dirac delta func-
tions on Cartesian grid points is feasible with appropriate interpolation schemes.
However, the overlap of grid points carrying redistributed charges and those involved
in the treatment Of geometric interface singularities leads to an accuracy reduction.
This interference becomes inevitable when a coarse mesh is pursued in an interface
treatment. Therefore, a Green’s function approach Of the singular charges becomes
attractive in our higher-order interface schemes for the linearized PB equation. Such
an approach decomposes the solution (2') into regular part (5 and singular part (17) [131].
The latter consists of a fundamental solution Of the Poisson equation with the singular

charge, (15*, and a harmonic function (1)0 [24]

where (13(r) is defined as

60‘) = (3-2)

28

Here ¢*(r) is the Green’s function

Nm

Mr) = Z ——3"——— (3.3)

z,:1€_[I‘ — I‘il

solving the Poisson equation with singular charges, while (1)0(r) is a harmonic function

in Q" satisfying the Laplace equation

V2600) = 0 in tr;
(3.4)

¢0(r) = —(j>*(r) on I‘.

The boundary value problem given by Eq. (3.4) is to be solved with designed order of
convergence over the irregular domain (2‘ with possibly geometric singularities. As a
solution to the PB equation, the singular part, (5(r), is solved subject to the following

jump conditions

[air = 0 and («my = -e‘V(¢* + a0) - nh‘. (3.5)

We work on the linearized version of Eq. (2.1), thus the equation for the correction

potential (13(r) = ¢(r) — C(r) is homogeneous
-V ' (€(I‘)V<P(l‘)) + K2(r)(<13(r) + PM) = 0, (3.6)

and can be rewritten as

—v - (6(r)va3(r)) + Hakim _—. —..2(r)a(r). (3.7)

29

This equation is solved with jump conditions

[PIP = 0 and I€¢~nlr = -I€<3nlr = é—VW‘ + (250) ' nlr- (3-8)

The interface problem given by Eq. (3.7) with interface jump conditions (3.8) and
appropriate far-field boundary conditions is to be solved with designed order of con-

vergence. The MIB method is generalized for this propose.

Solution to the boundary value problem

To obtain (1)0, we need to solve the boundary value problem given by Eq. (3.4). We use
the standard 7-point finite difference scheme on all the regular grid points. However,
for irregular grid points in Q”, the standard 7-point finite difference scheme requires
fictitious value(s) outside the interface. In this case, we could determine the fictitious
values by extrapolation involving the interface values. However, instead of using the
fictitious values outside the interface, we use the interface values at the intersecting
points of the mesh and the interface. Since the resulting finite difference weights Of
the 7-point scheme are not strictly symmetric, this treatment slightly reduces the
accuracy of the scheme at the grid points close to the interface. However, we have
tested that this local accuracy reduction does not affect the overall second order
convergence Of our scheme. Moreover, the resulting matrix is nearly symmetric and
diagonally dominant, and can be efficiently solved by a preconditioned bi-conjugate

gradient solver.

30

The evaluation of jump conditions

To solve Eq. (3.7) with the MIB method, we need to evaluate the jump conditions in
Eq. (3.8). In other words, after solving the boundary value problem (3.4), we need
to calculate the derivatives of d0 and d* in the normal direction at each intersecting
point of the mesh lines and the interface. According to Eq. (1.9), these normal
derivatives can be Obtained from the partial derivatives of d0 and d* in the 23-, y-,
and z-directions by means of the transformation matrix. We encounter the difficulties
that both d0 and d* are available only inside the interface F. In fact, the Green’s
function d"‘ is analytically available and its gradients can be computed to arbitrarily
high accuracy at each intersecting point. Therefore, the accuracy is Often limited by
approximating the derivatives of do. We seek to extend domain Q” by extrapolations
so that central difference schemes can be applied, since central difference schemes are
more accurate than the asymmetric ones for a given length of stencils. Moreover,
in many cases with geometric singularities and restricted geometric environment, a
domain extension becomes inevitable. Derivatives are computed on the enlarged
domain. Various techniques developed in our MIB are utilized to calculate these

derivatives. This procedure is discussed in the follows.

The extension Of d0 in the fictitious domain

Having solved the boundary value problem (3.4), d0 is available in (2’. To evaluate
partial derivatives of d0 at each intersecting point of the mesh lines and the interface,
we extend domain 9‘ by extrapolation. Fig. 3.2(a) gives three cases of a grid point
close to the interface and the schemes used to extend the domain. In Case I, there

are at least three grid points inside the interface in the direction of extension. We

31

Casel:

 

Caselk

Casein:

 

(a) (b)

Figure 3.2: d0 extension schemes (a) and illustration (b)

just use three d0 on grid points (2,3,4) and the value on the interface point B to
extrapolate the function value on grid point 5, and then use d0 on points (2, 3,4, B)
and the extrapolated function on grid point 5 to extrapolate the function value on
grid point 6. In Case II, there are only two grid points in between the interfaces
in the direction of extension. We use d0 on points (A, 3,4, B) to extrapolate d0 on
grid point 5, and then use d0 on points (A, 3,4, B) plus the extrapolated d0 on grid
point 5 to extrapolate d0 on grid point 6. In Case III, there is only one grid point in
between the interfaces in the direction of extension. We will use d0 on points (A, 4, B)
to extrapolate d0 on grid point 5. Due to the fact we only have three points here
and they are non-uniformly distributed, we will only extend d0 to one grid point to
maintain the accuracy of the extension.

There is a choice to be made when d0 on two extended grid points can be obtained
from extensions of different directions. It is found that the accuracy of the extension
is optimized if the following order of priorities is maintained: Case I for the first point,

Case II for the first point, Case III for the first point, Case I for the second point,

32

   

i i+4 i i+5

(a) (b)

Figure 3.3: Schemes to find partial derivatives: (a) regular cases; (b) an special case

and Case II for the second point. We have examined the accuracy Of this order by
using many test calculations with different geometric variations Of the interfaces. In
a smooth geometry with a sufficiently dense grid, Case I will be dominant and Case
III rarely occurs. Although our extensions involve schemes of different convergence
orders, these approaches are optimal or nearly optimal for geometric singularities,
which are extremely challenging. The overall scheme is of numerically second order
even for complex molecular surfaces of twenty four proteins.

After the extension, we obtain a fictitious domain outside the interface as illus-
trated in Fig. 3.2(b). The computations of partial derivatives a2, d3 and d2 are

carried out.

The calculation Of partial derivatives

To calculate the derivative of d0 in the normal direction, we need to compute three
partial derivatives, d2, d3 and d2, at each intersecting point of the mesh lines and
the interface. Since central schemes or nearly central schemes are preferred, these
derivatives are calculated by using d0 in the domain 90, on the interface P, and in
the extended fictitious domain. Fig. 3.3 gives an :L' — y cross section illustration of the

schemes for computing d200, j, k) and d300, j, k) at the intersecting point of the jth

33

mesh line and the interface, and near the grid point (2, j, k). In Fig. 3.3(a), the area
circled by the solid curve F is the domain (2‘, and the area between the solid curve
and the dashed curve is the extended domain. If we restrict ourselves to grid points
inside the domain 9‘, we can compute 2, by d0(2—2, j, k), d0(2—1,j, k) and d0(2',j, k)
but it will be impossible for us to calculate dg(2'0, j, k) directly. We therefore have to
consider d0 in the extended domain and on the interface. In this scheme, d2 can be
evaluated by using d0(2' — 1,j, k), d0(2',j, k), and d0(2' + 1,j, k). Similarly, dg(2o,j, k)
can be calculated by using d0(20,j— 1, k), d0(20, j, k) and d0(2'o,j+ 1, k) as denoted by
squares. However, d0(2'0, j — 1, k) and d0(2'0, j + 1, k) are not available and they are to
be computed by interpolations with (d0(2'—2,j— l, k), d0(2'— 1,j— 1, k), d0(2',j—1,k)),
and (d0(2'— 1,j+1,k),d0(2',j+1,k),d0(2+1,j +1,k)), respectively.

Very rarely but possibly, we may encounter more challenging situations as illus-
trated in Fig. 3.3(b). This typically occurs with geometric singularities. In this case,
d300, j, k) cannot be computed directly, nor can it be computed with d0(2'0, j — 1, k),
d0(2'0, j, k) and d0(2'0, j + 1,k) because d0(20, j — 1,k) cannot be computed even
with the extended fictitious domain. To overcome this difficulty, we propose an
alternative scheme. We first compute dg(2’, j, k) by using d0(2', j — 1, k), d0(2', j, k),
and d0(2',j + 1,k). We then calculate dg(2 + 1,j, k) by using d0(2' + 1,j — 1,k),
d0(2'+1,j,k), and 2% + 1,j + 1, Is). Finally, we use 2% + 2,j, k), 2% +2,j + 1, k),
and d0(2'+2, j+2, k) to compute dg(2+2, j, k). Eventually, d300, j, k) can be obtained

by interpolations with d3(2’,j, k), dg(2 + 1,j, k), and dg(2' + 2,j, k).

34

 

 

3.2.3 The computation Of electrostatic free energy of solva-
tion

The calculation of the electrostatic free energy of solvation plays a critical role in the
study of biomolecule in solvent. Such calculation is also one of the most important
applications of the PB equation. In the Green’s function formulation, the evaluation
of the electrostatic free energy of solvation is slightly different from that in the direct

solution of the PB equation

Nm
AG... = ;: grammar.) — Amie». (3.9)

i=1

where q(r,~) is the partial charges at r,- E Q', ¢diele and ¢homo are electrostatic
potentials in the dielectric and homogeneous environments, respectively. In MIBPB-II
and many other finite difference based methods, the singular charges are redistributed
to the neighboring grid points and the summation is computed over all of these
redistributed charges. Electrostatic potential in homogeneous environment, dhomo, is
obtained by the fast Fourier transform (FFT).

In the Green’s function formulation, the electronic potential in the homogeneous
environment, dhomo, can be identified with the Green’s function which solves the
Poisson equation with the singular charge source term, d*. Since d* is also a part of
the full solution as shown in Eq. (3.1), it does not contribute to the electrostatic free
energy of solvation. The latter, therefore, can be directly computed from the sum of

<5 and (:50
Nm
AG... = $2 quasar.) + 2%.». (3.10)

i=1

Note that d0 is defined in a smaller domain {2’ and its evaluation is relatively cheap.

35

3.3 Results and discussions

In this section, we examine the accuracy, validate the convergence and test the speed
Of the proposed MIBPB-III solver. For a comparison, we also applied some other
established methods, such as PBEQ [61], a representative finite difference PB solver
from CHARMM [18], and APBS [7, 52, 53], a recently developed multigrid finite
element and finite difference PB solver. The finite difference function of the APBS
is utilized in our calculations. The comparison is also given to our second generation
MIB based PB solver, the MIBPB-II, which provides a rigorous treatment of geo-
metric singularities. The specific MIBPB-II solver used in this comparison is a new
code. Extensive experiments are carried out over Kirkwood’s dielectric sphere [67],
molecular surfaces of one, two, and eighteen atoms, and a set of twenty four proteins
used in previous tests [34, 128]. Molecular surfaces are generated by using the MSMS
(v2.5.7) [102] at density 10 and with a probe radius of 1.4 A unless specified. In all
test cases, the dielectric constant is taken as 6+ = 80 and 6_ = 1. The latter is the
same as the dielectric constant in the vacuum, 60. For all protein structures hydrogen
atoms were added to obtain full all-atom models. Partial charges at atomic sites and
atomic van der Waals radii were taken from the CHARMM22 force field [64].

To compare the computational performance, we use three error measurements,
the surface maximum absolute error E1, the surface maximum relative error E2, and

the surface root mean square error E3

A

$(xiylz)—é(xtytz)
¢(x,y,z)

N- .
E3 = \/N1—‘ 27;:er [(152 — ¢2l2a

irr

 

E2 = maxr

 

 

 

36

where d and d are numerical and exact solutions, respectively. The maximum value
and summation are taken over all the irregular grid points near the interface I‘, whose

number is given by Ni”.

3.3.1 Units of the PB equation

The solution (I) of the PB equation given by Eq. (2.1) represents the electrostatic
potential, whose units could vary depending on which unit system will be adopted.
In this thesis, we use an intuitive way to define its units. We assume the model is under
molecular precision, therefore the application Of A as distance units will be an good
choice. Meanwhile, the charges carried by atoms are most conveniently measured by
ac, the charge carried by an electron. These two assumptions basically defines the
units of the values in our calculation. If we observe the right—hand side of Eq. (2.1),
the unit indicated will be ec/Il3 (The unit of the 6 function is 1/A3, since its integral
over the entire space gives constant 1). Consequently, we know that the units of the
potential will be given by ec/A if we use 13C and A as our basic units. However, the
most popular unit of potential as the output of many current bimolecular software is
given by kcal/mol - ec. From the calculation shown in Holst’s thesis, we need multiply
332.06364 (adjusted to 332.0716 to conform the calculation of CHARMM and APBS)
with our output d to give values of electrostatic potential with unit kcal/mol . cc.
All our reports in the following will use this unit for electrostatic potential and use

kcal/mol for solvation energy.

37

3.3.2 Validation

Kirkwood dielectric sphere

To validate the proposed MIBPB-III, we first consider Kirkwood’s dielectric sphere
of radius 2.0A. In our tests, the number of charges and their locations are allowed to
vary over a wide range. A detailed derivation of the model is given in the Appendix.
Case 1: A dielectric sphere with a centered charge.

We start our test with a single unit charge at the center of a sphere. For this case,

1

d* will be a function Of the radius 7‘, therefore it will have the constant value d* = F

on the interface, where b is the radius of the spherical biomolecule. This constant

_1_
be—

value of d* on the interface will lead to a constant d0 = — over the domain (2‘.
Since d0 is a constant, the interface jump condition (3.8) of Eq. (3.7) will depend
only on the normal derivative of d*, which is analytically given. Table 3.1 lists the
results of the potential near interface and the electrostatic free energies Of solvation
of MIBPB-III and two other established PB solvers, PBEQ and APBS. Results of
our previous MIBPB-II are also given to show the consistency and improvement of
the present MIBPB-III scheme. In this table, AG is the free energy of solvation with
exact value of -81.98 kcal/mol. All the solvation energies are reported in kcal/mol.
Meanwhile, all electrostatic potentials are reported in kcal/mol/eC unless specified.
At the mesh size of h = 0.5A, both MIBPB-II and MIBPB-III are able to provide
very accurate electrostatic free energy of solvation. In fact, their errors at h =
0.5A are smaller than those of APBS and PBEQ at the mesh size as small as h =
0.05A. These results indicate the impact and importance of rigorous mathematical

treatment of interface flux continuity conditions.

At the mesh size of h = 1.0A, MIBPB-II loses its accuracy in the electrostatic

38

free energy of solvation, because it does not have a Green’s function treatment of
charge singularity. In this case, the grid points that carry the information of interface
jump conditions are mixed with those that carry the distributed singular charges.
Consequently, MIBPB-II has a similar level Of errors as PBEQ and APBS, which do
not have an explicit mathematical treatment of interface flux jump conditions. In
contrast, the advantage of MIBPB-III is very obvious in this situation. Because of its
Green’s function treatment Of singular charge and its rigorous treatment of interface
flux jump conditions, the MIBPB-III is able to deliver high accuracy in the prediction
of the electrostatic free energy of solvation at h = 1.0A.

Surface electrostatic potentials are very important for protein-protein and protein-
ligand interaction, and thus, the electrostatic steering effects of signal transduction.
It was not widely recognized that traditional PB solvers produce very large relative
errors in the prediction of surface electrostatic potentials, as much as 70% at the mesh
size of h = 0.2A for a unit charge in a dielectric sphere. The earlier MIBPB-I and
MIBPB-II are able to deliver highly accurate surface electrostatic potentials for this
system, with less than 2% errors at h = 0.5A. However at a coarse mesh size of h =
1.0A, MIBPB-I and MIBPB-II cannot maintain their accuracy for the aforementioned
reason. Our MIBPB-III is designed to overcome this difficulty. Indeed, MIBPB-III
is able to provide extremely accurate surface electrostatic potentials at h = 1.0A as
shown in Table 3.1.

Case 2: A dielectric sphere with an Off-centered charge.

In Case 2, we allow the unit charge to locate at a distance a from the center
of the sphere. We set mesh size h to 1.0, 0.5, 0.2 and 0.1 A. Table 3.2 gives the
numerical errors of surface potentials calculated by PBEQ, MIBPB-II and MIBPB-

III. The unit of the surface potentials is kcal/mol/ec. From the table we can see that

39

both MIBPB-II and MIBPB-III can approximately achieve second order accuracy.
While MIBPB-III is more consistent in its convergence when the charge approaches
the boundary of the sphere. When both a and h are relatively small, MIBPB-III is
at least one order of magnitude more accurate than MIBPB-II. When a = 1.5 and
h 2 0.5A, MIBPB-II is abnormal because the system is under resolved. MIBPB-III
is still more accurate. In other cases, we see a pattern that MIBPB-III outperforms
MIBPB-II, and MIBPB~II outperforms PBEQ.

Table 3.3 provides a detailed analysis of MIBPB—III errors in Table 3.2. The
absolute errors are a factor of 332.0716 smaller in Table 3.3 because a different unit
for electrostatic potential (cc/A) is used. In this analysis, errors are computed for
different components, such as the flux, d0, d, and d. This analysis is important to our
understanding Of the advantage and disadvantage of MIBPB-III in more complicated
cases later. Note that Loo errors, which are measured over the whole domain, are
reported for d0 and d because they are well behaved functions. First, for a from
a = 0.2A to a = 0.6A, where the distance between the charge and the interface is
greater than 1.4A, it is seen that the flux, d0, d, and d can all approximately achieve
2nd order accuracy. For this distance range, most charges are located away from
the molecular surface. Second, for a from a = 0.8A to a = 1.2A, 2nd order accuracy
cannot be achieved for all these components and d. However, if we examine the errors
of do, d and d, it is seen that the errors at h =1.0A are very small, and very close
to those at h =0.5A, which results in the low order convergence from h =1.0A to
h =0.5A. This may be an advantage when we use the scheme on a coarse mesh like
h =1.0A. Third, when a is as large as 1.5A, where the distance between the charge
and the interface is about 0.5A, the order Of convergence of the flux, do, and d is

very low. However, even under this situation, the solution, d, is still of about 2nd

40

order accuracy. In short, although the accuracy of the overall solution is limited by
its components, it shows a reasonable trend Of second order convergence.
Case 3: A dielectric sphere with two Off-centered charges.

Our third case is a dielectric sphere with two off—centered charges located on the
x-axis and the z-axis, of a distance a from the origin. Table 3.4 gives the errors in the
electrostatic potential. The distribution of errors observed in this case is very similar
to that in Case 2. This study indicates that the number of charges does not have
much influence on the accuracy Of these numerical schemes. However, the distribution
of the charge(s) does have a major impact on the numerical accuracy. The problem
becomes more difficult when a is getting larger. In general, the proposed MIBPB-III
significantly outperforms other two methods.

Case 4: A dielectric sphere with six off—center charges.

Finally, we consider a sphere with six off-centered charges which are distributed
in two patterns. In the first pattern, denoted as Case 4(a), six charges are lo-
cated at (O.4,0.0,0.0), (0.0,0.8,0.0), (0.0,0.0,1.2), (0.0,0.0,-0.4), (-0.8,0.0,0.0), and (0.0,-
1 2,0.0). While in the second pattern, denoted as Case 4(b), six charges are located
at (0.2,0.2,0.2), (O.5,0.5,0.5), (O.8,0.8,0.8), (0.2,0.2,-0.2), (0.5,-0.5,0.5), and (08,-0.8,-
0.8). Some of the charges in Case 4(b) are closer to the interface. Table 3.5 gives the
error analysis for these two charge distributions. The proposed MIBPB-III achieves
second order accuracy. It is necessary to point out that at a coarse mesh size h = 1.0A,
when the charges are very close to the interface, i.e., within 1.0A, the computed po-
tential will not be as accurate as the case where charges are located far from the
interface. However, because the van der Vv’aals radii of the partially charged atoms
are normally larger than 1.2A, such small distance between the charge and the surface
will not occur in practical bimolecular simulations.

41

Table 3.5 also provides the electrostatic free energies of solvation. Case 4(a)
shows a better convergence in the solvation energies as the mesh is refined, because
its charges are located away from the interface. In contrast, for Case 4(b), larger
variations in the solvation energies are observed. However, these variations are smaller
than 2%, and are acceptable for most simulations. In the following subsection, we
test the performance of the proposed method for the combination of geometric and

charge singularities.

An analytically solvable model

Kirkwood’s solution in a dielectric sphere serves as a benchmark test for PB solvers.
However, it is also important to validate our method for its performance on com-
plex dielectric interfaces, particularly interfaces with geometric singularities. For this
purpose, we consider molecular surfaces, which are often non-smooth. One Of difli-
culties in testing numerical methods is the lack of exact solution for the PB equation
with complex geometry. Here, we construct an analytical solution for systems with
arbitrary singular charges and arbitrary dielectric interfaces. Consider the Poisson

equation

Nm
—V - c(r)\7d(r) + k(r) = 47r Zqi6(r — r,) (3.11)

i=1

with discontinuous {(1‘) across the interface P, which divides the whole domain into
two separated parts, the molecular subdomain Q" and the solvent subdomain 0+.
We set C— = 1 and 6+ = 80. In the singular source term, q, is 2th partial charge
located at r). Here, k = 3cosrrcosycos z. The problem in Eq. (3.11) is analytically

solvable and its exact solution is

42

N m

 

 

_ Qi

d (2:,y,z) = +coszcosycosz
2222-2 2 y.) .(._.,)2

¢+(:v,y.2) = 0 (3.12)

This is a valuable analytical test case and can be applied to PB solvers to check their
accuracy and convergence order. Obviously, it can be modified in many different ways.
In particular, it can be trivially modified for the linearized PB equation, namely,
adding a term, k2d, to both sides of Eq. (3.11) without changing the solution.

In Green’s function formulation, the solution can be decomposed into there com-
ponents, i.e. d = d* + do + d. Here, we can identify the solution to the Poisson

equation with the charge forcing term

N m.

.2 (12'
‘l’ 2222(2— 2'

2—1(—y y2)2 +(z-Zi)

 

 

As such, we have dO +d = cosxcosycosz in (2". Here, d0 is the solution of the
boundary value problem, Eq. (3.4), and d is the solution of Eq. (3.7). However,
due to the additional enforcing term of Eq. (3.12), the jump conditions Of the PB

equation are not given by (3.8), but are the following

[2] = 2* — 2‘ = — cosxcosycosz — 2* (3.13)
[625] = e+V2+ . n + e‘V(—2“ + 2* + 2°) - n (3-14)
[627,] = e+v2+ - m + e“V(-—2— + 2* + 2") - m (3.15)
[62.] = e+V2+ - 1+ e“V(—2‘ + 2* + 2°) -l. (3.16)

43

)T )T are the rows

where n = (P11.P12.P13 .m = (P21.P22.P23)T.and1= (P31,P32.P33

from the interface transform matrix p, given by Eq. (1.5). This benchmark test can
be used to validate detailed components of the numerical solution of the proposed
MIBPB-III, including its treatment of singular charges and its treatment of interface
singularities. In the following two subsections, we apply this test to our methods with
different interfaces, including the molecular surfaces Of polyatomic systems and the

molecular surfaces of proteins.

Molecular surfaces Of polyatomic systems

We first consider systems with one, two and eighteen atoms. For the monoatomic
system, a unit charge is located at the atomic center. The coordinates for the diatomic
system are (—2,0,0) and (2,0,0), with a unit charge at the center of each atom.
Atomic radii are set to 2.0 for these two systems. The molecular surface of the
diatomic system is constructed analytically. The atomic coordinates (1:, y, z) and van
der Waals radius (r) for the eighteen-atom system are given as (2:, y, z, 2"): (-2.0270
0.9540 —0.6510 1.7), (-1.6690 0.2340 0.6650 1.7), (—0.4530 -0.6870 0.4410 1.7), (0.7510
0.1480 -0.0400 1.7), (0.3930 0.8680 -1.3560 1.7), (-0.8230 1.7880 -1.1320 1.7), (-2.2840
0.2080 -l.4180 1.2), (-2.8880 1.6170 -0.4830 1.2), (-2.5270 -0.3680 0.9970 1.2), (-1.4260
0.9800 1.4350 1.2), (—0.1960 -1.1890 1.3850 1.2), (-0.7010 -1.4410 —O.3200 1.2), (1.0070
0.8940 0.7270 1.2), (1.6120 -0.5150 -0.2080 1.2), (0.1490 0.1210 -2.1260 1.2), (1.2510
1.4700 -1.6880 1.2), (—1.0810 2.2910 -2.0760 1.2), and (-0.5750 2.5430 -0.3710 1.2).
The partial charges of these atoms are artificially assigned to the atomic centers as
0.2 * (—1)', 2' = 1, 2...18. The molecular surfaces Of these systems are used in Eq.
(3.11) to test the proposed method. Note that the molecular surfaces of the diatomic

and 18-atom systems are very irregular. Both MIBPB-II and MIBPB-III have built

44

in rigorous interface treatments for these irregular features on the interface.

Table 3.6 shows both the numerical errors and convergence orders of MIBPB-II
and MIBPB-III. E2 cannot be measured for these cases because the solution vanishes
at some grid points. Second order accuracy is obtained by both methods. However,
MIBPB—III demonstrates a more consistent convergence pattern and shows a better
overall accuracy. Because of the large atomic radius of r = 2.0 in the monoatomic
and diatomic systems, we observe little interference effect of singular charges and
the interfaces. Therefore, two methods have the same level of accuracy at h = 1.0.
However, the difference in the performance of two methods at h = 1.0 is very dramatic
in the 18—atom system where some van der Waals radii are as small as h = 1.2. This

difference demonstrates the advantage of the proposed MIBPB-III technique.

Molecular surfaces Of twenty four proteins

It is well-known that molecular surface definition admits cusps and self-intersecting
singularities [102]. Therefore, it is important to verify that the proposed MIBPB-III
is of second order accuracy for the molecular surfaces of proteins. To this end, we
consider a set of 24 proteins used in our earlier tests [121, 128]. The molecular surfaces
of these proteins are used in our tests. For a comparison, our earlier MIBPB-II is also
employed. MIBPB—II has built in the rigorous treatment of interface singularities, but
without the Green’s function treatment for singular charges. Charges are distributed
to the neighboring grid points in MIBPB—II, as in many other finite difference based
methods.

Table 3.7 lists the numerical errors and convergence orders of MIBPB-II and
MIBPB-III for the molecular surfaces of 24 proteins. The proteins are listed in the

ascending sequence of the radius of gyration. Since numerical errors are collected

45

from the solution at irregular points around the interface, their magnitudes indicate
the reliability of PB solvers for the prediction of electrostatic surface potentials. The
most important feature in the numerical errors is that MIBPB-III is able to produce
much smaller errors at the coarse mesh size of h = 1.0A, which is crucial for the
proposed method being used in large scale computations. At h = 1.0A, MIBPB-II
is not accurate and reliable due to interference of charge singularities and interface
jump conditions. The accuracy of MIBPB-III at h = 1.021 is similar to or even better
than that of MIBPB-II at h = 0.50A. At fine mesh resolutions, i.e., h = 0.50 and
0.25A, MIBPB-III is about. an order of magnitude more accurate than MIBPB—II.
Moreover, MIBPB—III has a consistent order of convergence over two mesh refine-
ments. Whereas, MIBPB-II does not have second order convergence at the last mesh
refinement. This reduction in convergence order must be attributed to the presence
of the charge singularities, because the MIBPB-II was confirmed to have the second
order convergence for molecular surface singularities of these 24 proteins [121].

Fig. 3.4 depicts the maximum absolute errors (E1) at three mesh sizes. This
figure brings the above comparison in contrast with the errors of PBEQ. Obviously,
MIBPB-III errors are a few orders of magnitude smaller than those of PBEQ, which
are distributed between 100 and 101 for all meshes. The PBEQ fails to deliver accurate
solutions even if the mesh is further refined, because of the discontinuous nature Of

the solution at molecular surface interfaces.

3.3.3 Applications

Having validated the proposed MIBPB-III, we consider its application to the electro-

static free energies of solvation and electrostatic surface potentials.

46

 

 

 

 

 

8 A A AA 0 0 AA
A A A A. A
A A . O .0 g
O‘Oé¢6.ao'o.OOO-OO-ROOO A00 6
-2
10 r A
' A A :AAAAA
A AA AAAA A AAAAJ- A]
A t
-3
10 ‘
0 5 10 15 20 25

Figure 3.4: Comparison of absolute maximum errors (E1) of the surface potentials of
24 proteins, which are ordered with ascending radii Of gyration: 1ajj, 2pde, lvii, 2erl,
lcbn, 1bor, 1bbl, 1fca, luxc, lshl, 1mbg, 1ptq, 1vjw, 1fxd, 1r69, 1hpt, lbpi, 451C,
1a2s, 1frd, lsvr, 1neq, 1a63 and 1a7m. Errors of MIBPB-III (green), and PBEQ
(blue) at 24 proteins are respectively linked with dash-lines and solid-lines. Errors
of MIBPB-II (red) are unlinked. Here, El, 0, and A, are respectively for h = 1.0, 0.5
and 0.2521.

47

Electrostatic free energies of solvation Of proteins

Table 3.8 lists the details of electrostatic free energies of solvation of twenty four
proteins. For each method, results are reported at three mesh sizes, h = 1.0, 0.5 and
0.25A. Differences are computed with respect to the results Obtained at the finest mesh

2 0.25A, and their magnitudes can be regarded as an indication of convergence.
MIBPB-III and MIBPB-II have a similar level of accuracy and convergence at mesh
size h = 0.5A. However, MIBPB-II obviously produces larger differences at mesh
size h = 1.0A than MIBPB-III. In Table 3.8, the CPU time of two methods are
compared at h = 0.5A. At a given mesh size, MIBPB-II is faster than MIBPB-III due
to following two reasons. First, in solving the PB equation, due to the conversion of
the singular charges to the interface jump conditions, the matrix of MIBPB-III has a
relatively larger condition number and requires more CPU time. Second, to calculate
the electrostatic free energies, MIBPB-II solves an additional Poisson equation by
using the FFT scheme, which is very fast. Whereas, MIBPB-III has to solve an
additional boundary value problem.

Fig. 3.5 provides a comparison Of the MIBPB-III, PBEQ and APBS. First, there
is a good consistency in the solvation energies from three methods, as shown in Fig.
3.5(a). However, this consistency is illustrated at the scale of 3500 kcal/mol. In fact,
each method has a very different convergence property as shown in Fig. 3.5(b), where
the difference of solvation energies obtained at two different mesh sizes, h = 0.5 and
0.25, is plotted for each method. MIBPB-III demonstrates an excellent convergence
as the mesh is refined. Whereas, results of APBS and PBEQ show large variations
over the mesh refinement. As such, it is reasonable to regard the solvation energies

of MIBPB-III at mesh size h = 0.25A as converged ones and use them to assess the

48

 

 

—O--MIBPB-III (h=0.5)

+21,sz (h~0.25) 0 Wm
+PBEQ(h=O.25) .

kcal/mol

 

'70 +MIBPB—III
_80 —-<-—APBS

 

 

 

 

 

 

 

 

 

 
    
   

  
    
    

 

 

 

 

 

 

 

 

I.” PBEQ
"o 5 1o 15 20 25
(a) (b)
22 - - - 24c - . .
220_ +M|BPB-I|l (h=1.0)
o """ 200] —o—M|BPB-III (h=1.2)
—t>—Apes (n=o.25)
-20 180
40 160»
'g - -. g 14oL
. : C
> -so ~ . . 8 120-
E +MIBPB-III (h=0.5) ' 8 100-
'“ +M|BPB-III (h=1.0) so»
-o—MIBPB-III (11:12) _ 1
'100 +21sz (h~0.25) 60'
+Apas (II-0.5) 4o-
"120 +PBEQ(h=0.25) . 2°
+PBEQ(h=0.5) n
' ""o 5 1o 15 2o 25 "o 5 1o 15 20 25

c d)
Figure 3.5: Comparison of electrostatic free energies of solvation of 24 proteins listed
in the order of gyration radii. (a) Solvation energies; (b) Differences of solvation
energies between mesh sizes of h = 0.25 and h = 0.5A; (c) Differences of solvation
energies between MIBPB-III at h = 0.25A and the results in the legend; (d) CPU
time.

49

accuracy of other schemes. In Fig. 3.5(c), the differences of solvation energies between
those obtained by MIBPB—III at mesh size h = 0.25A and others are plotted. Indeed,
as the mesh is refined, the solvation energies of APBS and PBEQ converge toward
those of MIBPB-III Obtained at h = 0.25A. In fact, the magnitudes of variations of
MIBPB-III obtained at very coarse mesh sizes of h = 1.0 and 1.2A are in the same
range of those of APBS and PBEQ Obtained at h = 0.25A. Clearly, APBS and PBEQ
at the mesh size of h = 0.5A produce larger differences in solvation energies than does
the MIBPB-III at a very coarse mesh sizes of h = 1.0 and 1.2A. This analysis justifies
the CPU comparison of MIBPB-III at the mesh size of h = 1.0 and 1.2A with APBS
and PBEQ at mesh size of h = 0.25 for practical electrostatic calculations. As such,
we plot the CPU time used by APBS, which is well known for its fast speed, and
those used by MIBPB-III in Fig. 3.5(d). The CPU time is recorded on a Pentium-IV
PC with 2.8GHz CPU and 2G RAM. It is seen that at a similar level of accuracy,
MIBPB-III is about three times faster than APBS. However, it is necessary to point
out that at the same mesh size, APBS is much faster than MIBPB-III. Table 3.9

provides detailed data for the above comparison and discussion.

Electrostatic potentials

We consider the surface electrostatic potential of a heme-binding protein, Fe(II) cy-
tochrome C551 from the organism Pseudomonas acruginosa (PDB ID: 451C). The
heme group is removed and left with a cavity. A common computational domain,
[—17, 31] x {—37, 10] x [—9, 36], is used for all calculations so that computed electro-
static potentials can be compared on the same set of grids. The first row Of Table
3.10 lists the maximum and minimum of the electrostatic potential computed by us-

ing the MIBPB-III at mesh size h = 0.25A. These values provide an idea about the

50

variability Of electrostatic potentials on the mesh. The convergence of MIBPB-III is
studied by the differences of the electrostatic potentials computed by at mesh sizes
h = 0.5 and 1.0A, with respect to that computed at h = 0.25. There are relatively
small deviations in the potentials computed by MIBPB-III. In contrast, the differ-
ences of electrostatic potentials Of MIBPB-II Show a significant deviation. The same
feature is also Observed for PBEQ. One possible reason for these deviations is that
the grids that carry the distributed charges vary as the mesh is refined in these fi-
nite difference based methods. The differences of the electrostatic potentials with
respect to that obtained by MIBPB-III at h = 0.25 are listed in the last four rows
of Table 3.10. Dramatical variations are observed, indicating the mismatching of the
charge distributions in these methods. Note that the charge distribution is exact in
MIBPB-III, and the Green’s function solution for charges is also exact in MIBPB-III.
This finding may have substantial ramification for the reliability of the electrostatic
force computed by finite difference based PB solvers. A detailed study of this issue
is beyond the scope of this thesis.

The electrostatic compatibility of the protein and heme in the cavity region is
important to their binding aflinity. Fig. 3.6(a) illustrates surface electrostatic map-
ping of cytochrome C551 computed by MIBPB-III at h = 0.25A. The cavity region
is electronically positive except one site, which provides a favorable environment for
the heme. The convergence studies of three methods are given in Figs. 3.6(b), (c)
and (d), obtained by the difference of the potentials between mesh size h = 0.5 and
h = 0.25A for each method. It is seen that MIBPB-III shows an excellent conver-
gence at h = 0.5A, as there is little deviation in two potentials. MIBPB-II exhibits a
good convergence. Its deviations are evenly distributed on the interface, indicating its
errors are not created from geometric singularities. However, large deviations at the

51

 

(c) (d
Figure 3.6: Analysis of surface electrostatic potentials of cytochrome 0551. (a)
¢MIBPB_III(h=o.25); (b) ¢MIBPB-Ill(h=0.5) - ¢MIBPB—III(h=0.25)i (C) ¢MIBPB-II(h=O.5)
' ¢MIBPB-II(h=0.25)i (d) ¢PBEQ(h=0.5) ' ¢PBEQ(h=o.25)-

52

 

    

.51

( )
Figure 3.7: Differences of surface electrostatic potentials of cytochrome C551. (a)
¢MIBPB-III(h=1.0) - ¢MIBPB-III(h=o.25); (b) ¢MIBPB-II(h=o.5) - ¢MIBPB—III(h=0.25)i (C)
¢PBEQ(h=0.25) - ¢MIBPB-111(h=o.25); (d) ¢PBEQ(h=o.5) - ¢MIBPB—III(h=0.25)-

scale of :I:5 kcal/mol/ec are observed from the potentials computed by PBEQ, which
will have an impact in the electrostatic steering effect of the heme binding process.
The differences of surface electrostatic potentials obtained under various condi—
tions and that of MIBPB-III at h = 02521 are plotted in Fig. 3.7. It is seen that
electrostatic potential obtained by MIBPB-III at h = 1.0A shows considerable devi-
ations inside the cavity. The deviations obtained from MIBPB-II at h = 0.5A are
relatively small and evenly distributed. While the deviations obtained from PBEQ

are quite large, particularly for the coarse mesh h = 0.5A.

53

3.4 Conclusion

This chapter presents a novel Poisson-Boltzmann (PB) solver based on a Green’s
function formulation and the matched interface and boundary (MIB) method. Our
earlier MIB based PB solver, MIBPB-II, has built in an advanced mathematical treat-
ment of dielectric interfaces so that the subgrid information of the interface inside a
mesh is accounted and flux continuity conditions at the solvent-molecule interface is
rigorously enforced. Moreover, the MIBPB-II maintains its second order accuracy for
geometric singularities, including cusps and self-intersecting surfaces in the molecular
surfaces [71], without resorting to surface smoothing procedures. Consequently, the
MIBPB-II is more accurate at a coarse mesh than traditional PB solvers at a fine
mesh. Nevertheless, MIBPB-II suffers from an accuracy reduction if the mesh size is
about half of the van der Waals radius, because the interference of grid points that
carry the interface treatment and the grid points that carry the singular charges.
The present work effectively overcomes this difficulty by developing a Green’s func-
tion formalism for charge singularities. The essence is to split the solution into a
linear combination of three components in which the singular charge source term is
analytically resolved and leads to flux jump conditions at the interface for the PB
equation without the singular charge term. The MIB method is then employed to
solve the resulting PB equation and treat flux jump conditions from both the original
interface and the singular charge contributions on an equal footing. This new PB
solver, denoted as MIBPB-III, is extensively validated by benchmark tests, including
Kirkwood’s dielectric sphere [67], molecular surfaces Of one, two and eighteen atoms,
and molecular surfaces of twenty four proteins. Although molecular surfaces possess

geometric singularities, the designed second-order accuracy, i.e., numerical error re-

54

duces by a factor of four when the mesh is halved, is confirmed in the aforementioned
tests. It is found that the MIBPB-III provides some of the most accurate solutions
to the PB equation, to our knowledge. At a mesh as coarse as 1.2A, MIBPB-III is
able to deliver a similar level Of accuracy as other established PB solvers, includ-
ing PBEQ and APBS, at a fine mesh of 0.25A. As such, at the same accuracy, the
MIBPB-III is about three times faster than the state Of the art multigrid PB solver,
the APBS. However, at the same mesh, the MIBPB-III normally requires more CPU
time than PBEQ, APBS and MIBPB-II. Since the Green’s function solution to the
singular charges is exact, the present treatment provides correct peak values for the
electrostatic potential, which may be an advantage in the prediction of electrostatic

forces for molecular dynamics.

55

 

 

 

23. $2.... was- see was; was- was. an 8%- 3.8 an 8%.. Se

once was? 8.5- and $53 was- 8.8 £2.” 3%- 8.8 mm...” 22%. Se

83 names was- s; as...“ was- as: 3:. meme- 3.: 8s swam- one

See gees was- was one e.g- 85w 8.: 3%. Saw 8.: ES. one

83 ENE 3.5- 3% as meme- $8 3mm seem- 8.5 mass 3%. 2:

am cm 54 am am 04 am am 54 am 5 O4 c
32592 22:55 8%. 3mm

 

 

 

 

 

S $.de 221$. oEomEmC e E mmneno E: e. “8 Echo EEQBOQ 835m 98 mowmaoco deb Beaumoboflo mo :OmfimQEoU ”H.” Each.

56

 

 

 

 

 

 

 

 

57

 

.2: 3333 mi 2133 :3 2:33 3
z: 8&8; :3 S+E§ 23. 853.3 3
33 21mm: «3. 83:3 9;. 8+3: 3
85:3 8+3; 8+3? 3 2
m3 833 3m 85%; 3.3 8533 3
N: 833.“ mm: 8+mm3 :3 8+3? 3
83 8+me 8.7 8.323 3? 8+3: 3
8534 8553 83:3 3 3
33 8.33 :3 833 m3 8+3? 3
2: 838.3 33 33mg.” :3 2:33 3
m3 833. :3 8+3? 33- 8+3: 3
Sums: 2:53 SE33 3 3
Q3 8333 m3 83m.” :3 8+3? 3
m3 Sums; 35 Rims.“ 33 S+mm§ 3
N3 533 33 8+3? o3- 8+3? 3
8.33 8+3: 85$.” 3 3
«3 33.33 :3 8.33 33 8+3; 3
a: 833 w? 8+9: o3 8+3: 3
m3 833. m3 aims.” 83 2:33 3
833 S+E§ 2133 3 3
v: 833 m3. 833 ”3 85.8.3 3
o: 833 m3 832 o3 8+3: 3
m3 833. 2: 2133 33 8+3: 3
533 83:2 8+m£3 S 3
x: moms.“ :3 8.33.3 m3 8+3; 3
m: Eats 83 833 33 853.3 3
$3 3333 33 8333 33 5+3; 3
333 8&8; SE83 3 3
£5 E 820 5 5E0 5 a 3
3-332 2-95:2 3mm

 

 

 

 

.Am manV 221% 2.3.536 w E v.33? wmumusmufio 5w “8 fldgauom 0353588 9:53 as $0.29 3258: Z ”Nd wEmH.

 

 

 

 

 

 

 

 

 

nwA moflAdA wnA 3&th 3A mdumdmA mmd- Adm—3A Ad
AdA Mdumnwd ddA mdnmNmA 3A MdAmAmd mdN AdAm—mdA Nd
wmd NdfiwvA mdN- MdAmdNAV wvd NdfiwmA mdd Adumwwd nd
NdLmAmdA mdAmmoA NdAmmoA dd+mAmA dA mA
NmN VQQEA NAN mdflddd 3N «dfldmA dmA modawmd Ad
NNA voided mdA voAmNdd EA wofidms mwA NdnmmAN Nd
wvd MdumtlAAV dAd MdummmA wmd 8&de wAA AdflbAA md
MdumAwd 8&va moaned AdfimdN dA NA
3N mdumnvd wNN mdLmAmNN me 3&de NNA mourmNdN Ad
2: voflmvN mhA voflmoA de VdnmAddA NwA moflNdN Nd
Nnd monvA Add- vofiovh cwd 8&va mmA Nodmud md
mdumdmN gamma. mdfimvN Ndéwmd dA dA
mmN mduvaA EA mdAmdAA de ddAmes omA woflvod Ad
de mdfimuh mmA mdflwnd wNN mdAmNod mwA mofldAN Nd
Ndd voflwmd dmA demNNN mAA weaved wvA NdlmANAA md
3&va nwoumAANnm 3&de NdflNAd dA wd
3A mdAmddA dmA mdlmddA EN ddnmdmA VAN onANA Ad
AdA mdfiwdd ANA mdflmAd AVN mdlmmdA mwA wofivmd Nd
dNN voAmmmA dN.N wofldmA dmA mdummmd NdA mdAmNdN md
voAmddd voAmAmd vdAmAmN mdlmddd dA dd
EA ddlmde EA ddflmmd EN SQNAN AKA mdfith Ad
EA mdfiwhN NhA mdumduN AmN mod—3A wmA mdflhwd Nd
mmN onNmA mmN wed—3A hmA mdumva mnA fidummmd md
vdAmNnd Audkmtd mdLmANm...V mofiuwA dA vd
va ddflgd va mod—3d va mdAmANd AKA ddAmNNd Ad
AKA mdumdAN EA mdeAN ddN mdAmAdd mmA mdflnoA Nd
3N Adm—SA mvN wofiddA dAN bdummAN 5A mdumddav md
vo‘mAwd vdfiAwd ddummdd woubwoA dA Nd
EEO Am EEO 8A $20 8A SEQ Am a 6
She Av Hobo mv hobo 9% Hobo USE

 

 

 

 

 

.AN mamov 823m 2502va a E mméco vmgvuzmofio SW pod £33303 0538.583 83:5 do £94.25 quGOQEoO ”m.” 29¢.

58

Table 3.4: Numerical errors of surface electrostatic potentials for two off-centered
charges in a dielectric sphere (Case 3).

 

 

 

PBEQ MIBPB-II MIBPB—III
a h E1 Order E1 Order E1 Order
0.2 1.0 5.91E+01 1.82E+01 3.72E-01

0.5 4.14E+01 0.51 1.02E+01 0.84 5.44E—02 2.77

0.2 1.81E+01 0.90 7.80E—01 2.81 1.17E—02 1.68

0.1 9.14E+00 0.99 1.00E—01 2.96 3.26E—03 1.84
0.4 1.0 6.67E+01 2.07E+01 4.32E—01

0.5 4.94E+01 0.43 1.16E+01 0.83 6.58E—02 2.71

0.2 2.16E+01 0.90 1.01E+00 2.66 1.50E—02 1.61

0.1 1.09E+01 0.99 1.10E—01 3.20 4.39E—03 1.77
0.6 1.0 7.17E+01 2.05E+01 4.22E—01

0.5 5.73E+01 0.32 1.27E+01 0.69 6.12E—02 2.78

0.2 2.60E+01 0.86 1.38E+00 2.42 1.65E—02 1.43

0.1 1.33E+01 0.97 1.80E-01 2.94 5.17E—03 1.67
0.8 1.0 7.24E+01 1.61E+01 4.88E—01

0.5 6.41E+01 0.18 1.04E+01 0.63 1.69E—01 1.53

0.2 3.16E+01 0.77 2.12E+00 1.74 2.13E—02 2.26

0.1 1.69E+01 0.90 3.00E—01 2.82 5.36E—03 1.99
1.0 1.0 6.74E+01 5.01E+01 1.59E+00

0.5 6.79E+01 -0.01 1.91E+01 1.39 5.96E-01 1.41

0.2 3.92E+01 0.60 3.58E+00 1.83 8.57E—02 2.12

0.1 2.23E+01 0.81 5.80E—01 2.63 1.46E—02 2.55
1.2 1.0 5.81E+01 2.04E+01 3.86E+00

0.5 6.70E+01 -0.21 3.17E+01 -0.63 1.61E+00 1.26

0.2 4.87E+01 0.35 5.80E+00 1.85 3.00E—01 1.83

0.1 3.10E+01 0.65 1.34E+00 2.11 5.94E-02 2.34
1.5 1.0 3.72E+01 7.52E+00 1.28E+01

0.5 5.00E+01 -0.43 4.28E+01 -2.51 5.40E+00 1.25

0.2 5.95E+01 -0.19 2.48E+01 0.60 1.98E+00 1.09

0.1 5.15E+01 0.21 6.37E+00 1.96 5.39E—01 1.88

 

 

 

 

 

 

 

 

 

 

Table 3.5: Numerical errors of surface electrostatic potentials for six charges in a
dielectric sphere (Case 4).

 

Example h E1 Order AG
(3.) 1.0 3.25 —3011.12
0.5 1.72 0.918 -2995.41
0.2 3.09e-1 1.87 -2990.20
0.1 5.91e-2 2.39 -2989.52
(b) 1.0 20.2 -3177.64
0.5 4.68 2.11 -3159.33
0.2 5.84e—1 2.27 -3122.70
0.1 9.46e—2 2.63 -3123.94

 

 

 

 

 

 

 

59

 

 

 

 

 

 

 

and 8N a; m: wow «3 520
.3de 3mm: norms; 8&2 Sums.“ 8de gamma Samoa E
as 3m RN mm.“ 3.. EN 35
8333 games 8%? Norma; Edge Ease Semen same” a mass M:
mm.“ :3 EN SN 8.” 8a 520
8&2...” 3am“: gamma Baas SEE.“ media 853 momma... E
was 9..” as an... 2.” ”3 520
Bass 3&8.” summed saw: 333 833 833 Saga 2 was:
8.“ new as o: a.“ 8.“ $20
8E3 3%: SEE Sumo? morass Sums.» 883 saw: E
mum . a: 9: new and a: $20
833 Edam...” modem.“ 8.sz worms.” normed meats News.” a 8% H
33 mg no 3 mm; was 3 .3 mama: magma
mm 5

 

 

 

 

manta 52532: SE.» A29 .5 $58 H8 Adrmmmda was BimEE do 853968 EB meotm dd 03a?

60

 

 

 

a: 8.888 :8 8888 8+5? 88 8.88.8. 88 8888 88:3 :83
8.8 88:88 8.8. 8888 8+m88 88 8.83 :8 8888 888.8 82
88 88888 83. 8888 8+8? 888 8.8.8.8 8.8 8888 88.88 as:
88.8 8888.8 88 88:3 8+888 88 888.8 88 8.888 8888 .92
88 88888 88. 8.8.88 8+8? 88 8888.8 88.8 8888 8888 E:
888 88:88 888 8888 8888 888 8888 88 8.888 88:88 .83
83 888.8. 8:. 88.88 8+m88 88 8888 888 8.88.8 88.8.8 28
888 88888 88 88:88 88:88 88 8888 88 88:8 8888 8.:
88 8888 $8 8888 21888 888 8888 88 888.8 8888 88
83 8888 88 8888 8+m88 88 8.88.8. :8 8888 8888 8:
8.8 88888 88 8888 8+m38 88 8888.8 8.8 8888 8888 88:
8.8 8888 888 888.8 8588.8 :8 8888 88 8888 88888 3.5
88 88888 8.8 8888 88.83 :8 888.8. 88 888.8 8888 33
82 88:88 88 8888 8588 88 8.888 88 8888 8.88.8 8888
32 8.88 8.2. 88888 8+8: 8.8 88888 888 888.8 88.83 :3
888 88.88 88 8888 85.888 :8 8888 88 88888 8888 an:
88 8888.8 88.8. 88.888 8+m888 B8 8888 88 888.8 888.8 8:
88.8 88888 88 888.8 88E: 88.8 8888 88 8888 888.8 38
8.8 8.83 8.8 88888 8+m88 88 8888 888 8888 8.888 .88
88 8.8.3 :8 88.88 888.8 88 888.... 888 8.8.8.8 8888 32
88 8.888 88.8. 88888 8+8? 2.8 8888 8: 88:8 8888 88
88.8 8.888 88 88:88 88:8 38 8888.8. 888 8888 8888 83
888 8888 88 888.8 888.8 8.8 88.8 8.8 8.8888 8888 288
82 8.888 88 88888 8888 88 8888 88 8.888 888.8 .83.
85 888 88.6 88 3 EEO 888 EEO m8 3 53on
2-382 5-88:2

 

 

 

85305

Sod .3283 do 888258 8:529: 23 58> AAAdV .Um 9:38 8% EWmnAmEZ 8:8 E-mnAmAE do 8828385880 28 88.8 AM RE. "82m GENE

61

 

 

 

88 88 8.3.- 8.888- 2- 8.88- 8.88- 888- «8:8. 88 8828- 8.88- 88-2
88 88 8.88- 8888- 88- 8.888- 88888- 8.8- 8828- 8.8 8.888- 8888- 82
888 888 8.8. 88:- 88- 8.8:- 8885- 8.8- 8.8:- 88 E8:- 885- 8E
888 888 888- 8.8:- 88- 888:- 8.2:- 888- 8.8.:- 88 8.8:- 8.2:- :2
SH 88 8.8- 8888- 3- 8.888- :88- 888- .8888- 88 8.888- 3888- E:
E 88 8.8- 8.82- 3 88::- 8888- 8.8- 8882- : 8.28- 8888- 82
8 88 888- 8.888- 8- 8.888- 8.888- 3 8888- 8.8 8888- 3.88- 28
S. 888 8.8- 8888- 3 8.888- 8.88- 8.:- 8887 88 8.888- 8.88. fig:
8 8: 8.8- 8.88- 8.8 8888- 888- 38- 888- 88 8.28- 8.38- 35
8 8: 888- 8.2:- : 8888- 8.88- 888- 8.2:- 88 8888- 8.88- 8:.
8 88 8.8- 8.888- 88 8.888- 8888- 8.8- 8.888- 88 8.888- 8.888- E:
8 82 8.8- 5888- 88 88888- 8888- 88- 8.888- 88 88888- 88888- a?
E- 8 E:- 8888- 88 38- .888- 88- 88- 88 8888- 8888- 33
8 88 888- 8.888- 3- 8.888- 8888- 8.2- 8.888- 3- 8.388- 8.888. 88::
8 8 888- 88:- 88- 888- 888- 8.8 883. 3 888- 888- 3.0.8
8 8 888- 8.8:- 88- 38:- 8882- 8.3- 8.8:- 88 888:- 8882- an;
S 88 8.88- 8.88- 88- 8.888- 8.888. 8.8- 3:88- I 8.8:- 8.888- 8.:
8 8 8.8. 888- 88- 888- 888- 88- 888- +8 888- 888- 35
8 8: 888- 8.888. 8.8- 8.88- 8.88- 88 :8- 88- 8.88- 8.88- 85
8 8. 8.8- 8.88- 8.8.. 8.88- 8.88- 88- 8.88- 88 8.88.. 8.88- 52
88 8 8.3.- 888- 8.8- 888- 888- 888- 888- 88 888- 888- E8
8 8 88- 8.88- :- 888- 8.88- 8.8 8.88- 88 8.88- 8.88- a:
8 88 88- 888- 88 8.88. 8888- 88 888- 8.8 2:8- 8.88- 888
88 8 888- 8.8:- 88- 8.3:- 888:- 3 H88:- 88- 888:- 888:- .53
= E .85 8.8 .85 88 888 .85 8.8 .85 88 888 £828
:5 8888882 58882

 

 

 

 

85580.5 58 3:95 we 503838 We 82385 ooh 288880585 "w.” @388

62

 

 

 

Sm E 5: 38m- 32m- Sam. «.85- Scam. EEN- 32m- 32m- 52
a: on E Sag- 33m- @83- Sovm- 3mg.- m.m€m- 35“- flumm- 83
3 mm mm 5.8m:- ofit- 25:. 3E. 5::- 352- EN:- SE- 83
om on 3 3w:- odET SE- 33- EE- hmvz- 38:. EE- :3
mm cm on 3am- 3%? 92mm- ©.Ewm- 38m- Emma- flamm- mSwm- E:
2: mm S 9&2- 3mm:- 322- @82- 33? 3mm:- wfia- 022. mg
on mm mm 3on- OQS- came- came- vde- 3%:- SNS- 32:. £3
8 a cm .932- 35- QEE- 3:2- vfifl- 3::- 232. @82- E:
on S cm 3%- mdmw- 25. 3%- 3;. m.mmw- $5- mix- pg:
3 2 mm 6%:- m.mo:- 92-:- wde- 35:- m2:- «.82- 0%,:- 8:
S E 2 53m- 33m- 08%. vfimm- 93%- Sam. 25mm- 98mm- 3:
g 2 S 287 3mm? @82- mgfi- 3mm? 3va- osmfi- EMS- is
S 2 M: ema- Emw- has- 3%. cdww- odww- 3%. 3%. SE
8 S cm 387 m8»:- HSS- $5- 532- $3:- mfié- H.32- was:
S a 2 «new- Nak- 38- $3- SE- SE- Qfi- 2E- 32
:- 3 S 35:. 32? Em:- mSS- SE:- sz- mt“:- sz- 8a:
B 2 2 rag-S- v.82- wfi-S- 3::- Eofi- 3:2- odaz- 252. .8:
an 2 2 $2:- mdoS- ~32- 3%- 32:- 33- m.©wm- memo.- BE
3 E 2 3%- 3%. 3%- mdow- 3mm- :3. 3mm- Emm- Hon:
3 S B QR”- wém- 93.“- odom- ”.mmm- m8”- 3%. w.mom- 52
:- m 2 3mm- ;8- 3%. 38- 33- 38- 3g- £3- :8
3 m 3 meme. 38- 3mm- mdoa- wim- mgw- $8- 38- a:
9- 2 2 $3M- 3mm- 3mm- 3%- :a. 32w- 2%- mdmw- £3
9- w 2 .69.:- 337 So:- 337 532- 33. 3”:- mst. .3:
3o 2 3 3 30 0o 30 2 3 0o mg Q
mmnE 5.95% SE mmn? E-mmmzz

 

 

 

Amwzogmv 25-H- Dav

 

CoE>wuxv Epmcm :oniom

 

 

9:805 58 5:35 we Sosa-{Lem *0 mwflwgmzw 8d osgmobuflm ”a.” @358

63

Table 3.10: Electrostatic potential values and differences (kcal/mol/ec) for cy-
tochrome C551.

 

 

 

 

 

Methods I Maximum value Minimum value
éMIBPBJIILhZOQQ 8230.13 -5515.96
CDMIBPB-uuhzos) - ‘Z’MIBPB-Ill(h=0.25) 34-54 -38-11
QbMIBPB-Ill(h=1.0) ' ¢MIBPBJII(h=0.25) 150-64 -45-68
¢MIBPB-Il(h=0.5) ‘ ¢MIBPB—Il(h=o.25) 980-89 ~968-10
¢PBEQ(h=o.5) ' ¢Peeggh=o2m 984-64 -967.67
¢MIBPB-II(h=0.5) - ¢MIBPB-Ill(h=0.25) 2490-40 4242-73
¢MIBPB-II(h——-O.25) - ¢MIBPB-III(h—_—0.25) 4526-79 43428-70
¢PBEQ(h=0.5) - ¢MIBPB-III(h=0.25) 2491-59 4242-03
¢PBEQ(h=0.25) ' ¢MIBPB-III(h=0.25) 4525-75 45428.72

 

 

64

Chapter 4

MIB Based Molecular Dynamics

Method

4.1 Introduction

The applications of the PB equation include electrostatic solvation free energies
[37, 70, 107], pKa values [43, 92, 117], and electrostatic analysis of interaction sur—
faces. In the past two decades, there has been continuous effort in developing PB
equation based electrostatic force calculations for molecular dynamics (MD) [45, 84—
86, 97, 105]. Theoretical foundation for such calculations was established by Gilson
et al [45] based on the free energy functional derived by Sharp and Honig [105]. To
stabilize the electrostatic force calculation involving the solvent-solute interfaces, Im
et al. presented a smoothed dielectric boundary approach [61]. Using a smooth
surface approach, Lu and Lou demonstrated that PB equation based methods can
be successfully and consistently applied to stable and efficient dynamics simulations
of unconstrained biomolecules [85]. Remarkably, Prabhu et al. reported a smooth-

permittivity finite difference PB method that is up to eight times faster than the

65

explicit water simulations [97].

Both Poisson-Boltzmann (PB) and GB models utilize a solvent-molecule interface
to separate the continuum domain from that of the biomolecule. Various biomolecular
surface models, including the van der Waals surface, the solvent accessible surface
[71], the molecular surface (MS) [99] and minimal molecular surface [11] have been
employed in the PB methods. At a coarse level of approximation, these surfaces
produce similar results for many PB applications. Thus, the selection of a surface is
often dictated by the numerical convenience. In particular, sharp interfaces, such as
molecular surfaces are often avoided in PB based MD simulations because of possible
numerical instability [85]. Smoothing techniques were used to avoid this problem
by producing smoothed surfaces [48, 61, 85, 125]. However, Swanson et al. have
shown that some atomic centered dielectric functions admit unphysical high dielectric
interstitial regions in implicit solvent models [108] and may lead to 40% overestimation
of solvation free energy [108, 109, 115]. Consequently, classic molecular surfaces and
new minimal molecular surfaces [11] are preferred choices for the PB theory. Due to
numerical instability, non molecular surface based solvation force calculation has been
reported. The accuracy, stability and reliability of existing PB based MD methods
have not been systamically validated in the lierature.

The objective of the present work is to develop MIB based stable and accurate
algorithms for electrostatic force calculation and molecular dynamics based on the
PB equation and molecular surfaces. The MIBPB technique is utilized to handle
molecular surfaces. There are a number of new features in the present algorithm.
The present force calculation is based on the molecular surface [99] or the minimal
molecular surface [11], without resorting to any surface preprocessing, such as the ge-

ometric smoothening and / or the dielectric smoothening. These biomolecular surfaces

66

could avoid unphysical high dielectric interstitial regions [108]. However, technically,
a few modifications to the existing force evaluation procedures are required. First,
force density expressions, particularly the dielectric boundary force density, need to
be modified. Second, the resulting force density involves the evaluation of limiting
values of the first order derivatives of the electrostatic potential from both inside
and outside the interface. Third, the molecular surface interface contains reentrant
surfaces which are detached from individual atomic spheres. Therefore, forces from
reentrant surfaces need to be distributed to nearby atoms. Fourth, accurate surface
integrations are required for the evaluation of dielectric boundary forces and ionic
forces. These evaluations are to be carried out along the complex interface and in the
Cartesian grid. Finally, a PB solver that is able to handle sharp interface without
numerical instability is needed. This is done with our MIBPB-II method.

The rest of this chapter is organized as the follows. Theory and algorithm are
presented in Section 4.2. The derivation of electrostatic force density functions is
based on a free energy functional. Modification of the density functions is discussed
with respect to sharp interfaces and molecular surfaces. A number of numerical
algorithms, including the evaluation of the force density as limiting values at the
interface, accurate surface integrations in Cartesian grids, and the distribution of
dielectric boundary forces to the individual atoms, are developed. Validation and nu-
merical experiments are described in Section 4.3. The method of accurate integration
in Cartesian grid is tested with the area of molecular surfaces. The proposed force
algorithms are tested over many polyatomic systems. Comparison is given to thin
the CHARMM PB based MD method, the PBEQ [61]. Section 4.4 is devoted to the
molecular dynamics applications. The present electrostatic forces are implemented in

association with the AMBER force fields, particularly, in place of the AMBER PB

67

module [85, 86] for molecular dynamics simulations. Two small biomolecules are used
in our demonstration. This chapter ends with concluding remarks summarizing the

results.

4.2 Theory and algorithm

4.2.1 Free energy functional

Consider an open domain 0 E R3. Let I‘ be the interface which divides 9 into two
disjoint open subdomains, Q = Q" U 0+. Here 9‘ is the biomolecular domain, and
9+ is the solvent domain. The starting point for the implicit solvent theory is a given
set of atomic coordinates for a molecule in R3. The charge density of the molecule
pm(r) can therefore be obtained by assigning charge to each atom according to its
type and environment
Nm
Pm = 47f ZQi6(r — 13)-
i=1
where Q,- is a partial charge, Nm of the number of charged atoms. The positions of
the charge r, is determined from experiments, such as X-ray crystallography, nuclear
magnetic resonance and computations. The permittivity €(r)
e— r E Q“ U P;
6(r) = (4.1)
6+ r 6 0+.
is a function of position r E R3, and its value in the molecular domain (9') differs

from that in the solvent domain (0+). The charge density of the solvent p3(q‘)(r))

is a function of the electrostatic potential p3. In case of a multi-species system, a

68

commonly used ionic density is the Boltzmann distribution and is given by
N ab/kT
p3 = AZqJ-cje q] , (4.2)
1 =1
where N is the number of ionic species, cj’s are the bulk concentration of each ionic
species, and qj’s are charges of each ionic species. Here )- is a characteristic function
indicating the domain of the solvent

0 r e Q— U P;
Mr) 2 (4.3)

1rEQ+.

The free energy of the electrostatic system can be expressed as an integral [105]

Gelec : fgelec(xi ya Z, $3 @137 (by: ¢Z)dQ: (44)
Q
where g918C is the free energy functional or energy density

6
9818C = P1an _ '2‘(V9b)2 + [33- (4-5)
Here [3,-(9‘)) = fp3(qb)dqb and [35(0) 2 0. If p3(¢>) is given by Eq. (4.2), then we have
N
[3,. = —kTZ[c,(e-q2:¢/kT —1)])-. (4.6)
i.

To minimize the free energy, the energy functional gelec(1:,y,z,¢,qu,¢y,<z5z) must

satisfy the Euler-Lagrange equation

agelec a agelec a agelec a agelec

——-—— —— —— —— _ 4.
(9a (9:0 00;.- + 8;, (9a,, + 82 (9qu ’ ( 7)

69

which gives rise to the PB equation,
—V - €V¢(r) 2 pm + p3. (4.8)

Obviously, there are flexibilities in the selection of p3. In particular, the present choice
does not account for the effect of ionic correlation, and ion-macromolecule interaction

near the interface.

4.2.2 Electrostatic forces

The electrostatic force Felec can be expressed as the gradient of the free energy density

elec

9

Felec = —/Vgelecdf2. (4.9)
Q

In this treatment, we assume that the system is at equilibrium so that the temperature
T, the ion concentration c,- and their charge q,- are position independent quantities.

From Eqs. (4.5) and (4.6), we have

N
Felec .—. —/ ($me—$(V¢)2V£—kTZ[Q(e“qi¢/kT—1)]VA (4.10)
9 i

N
+ pchD — cw) - VV¢ + )- Z[q,c,e-qi¢/’~‘Tv¢] d9.

1

The first term on the right hand side of Eq. (4.10.) can be integrated by parts

/ ¢medQ= - / pmVon= /(me)d§2, (4.11)
Q Q

70

where the product (ppm vanishes at infinity. Here E = —V¢ is the electrostatic field

vector. Using the product rule on the fifth term of Eq. (4.10), we have

fave) - (vvamn = / v . (New-mo— /(V¢)V - (chb)dQ. (412)
Q Q Q

The first term on the right vanishes as a surface integral at infinity. The resulting

force expression becomes

N

F8160 = — / me— $E2Ve—kT2 :c,(e-qz'¢/kT—1)VA (4.13)
D .
2

N
+ pm + V - (€VQ5) + AZqicie’qi‘p/kT (ng) d9.
i

Recognizing that the terms in the square bracket vanish because of the PB equation

(4.8), one arrives at the force expression

N
Felec = — / me— $E2Ve — kTZc,(e—qi¢/kT —1)V)\) an (4.14)

1-

Q

— felecdfl
Q

where the force density is [45]

N
fflee = me — $E2V€ — kTZc,(e_qi¢/kT — 1)V)\. (4.15)

1-

Here we have used the sign convention of Gilson et a1. [45]. The total force can be

broken into three contributions Felee = FRF + FDI3 + FIB. The first term is the

71

reaction field force

FRF = —- / medo (4.16)
9

Nm
= 47-; [a Q-V¢6(r — r.)

Nm

2 47f 2 Qt (qu) lr=ri ~

2:]

FRF

In practical computations, corresponding to the solvation is evaluated as the

difference of the solutions in solvent and in vacuum

FRF

Nm
471' Z QillVQbsllrzri — (Vfibvllrzril (4-17)

i=1

Nm
= :Fa
i=1

where F,RF is the reaction field force at 2th atom. The second term in Eq. (4.14) is

the dielectric boundary force

FDB = / éEZV-sdo. (4.18)
0

The last term in Eq. (4.14) is the ionic boundary force

N
FIB = / krzc,(e—qz¢/kT—1)wdo. (4.19)
9 2'

72

4.2.3 Force expressions involving sharp interfaces

For sharp interfaces, the permittivity values 6(1‘) at molecular domain 52‘ and the
solvent domain 9+ are different

€_ r E Q— U P;

6(r) = (4.20)

6+ r 6 9+.
In most calculations, 6’ and 6+ are set to 1 and 80, respectively. This jump in
permittivity gives rise to discontinuous coefficients in the PB equation (4.8). To
ensure the uniqueness of the solution, the PB equation is to be solved with the

following continuity conditions at the interface

[<25] = ¢+(r)-¢_(r)- (4-21)

[cqbn] 6+(r)V¢+ (r) - n — e—(r)V¢_(r) - n, (4.22)

where n = (7133,7131, nz) is the outer normal direction of the interface I‘. Here [qfi] = 0
for the PB equation and [6%] are in general non-zero. The interface normal can be
computed from given interfaces, such as the molecular surface of a given protein.
For the dielectric boundary force, although the V6 can be treated as a delta
function evaluated at the interface, the electrostatic field vector E = —V¢ is undefined
at the interface because of the flux continuity condition Eq. (4.22). Therefore, the
dielectric boundary force and ionic boundary force provided by Gilson et a1 cannot be
directly used and have to be modified. To this end, we express the dielectric function

as

e = e+H+(n — R+(r)) + e—H"(n — R‘ (1‘)), (4.23)

where H + is the Heaviside step function, H ‘ is the reversed Heaviside step function,

n is the normal vector that passes through an interface point R(P), and

Ria‘) = lim R0“) :1: efi, (4.24)

e——>0

where 6 is the unit normal vector. Since VHi(n - Ri(I‘)) = ifi6(n — Ri(f‘)), the

dielectric boundary force is

FDB = / / $1212 [#66- — R+(r)) — 575(7- — R-(r))] fidndS. (4.25)
P 11

By carrying out the integration over 71, we have

FDB = / fDBdS . (4.26)
I‘

2 fr %(€+|E+)2) — e—IE‘|2)dS, (4.27)

where dS = fidS is pointed at the surface normal direction 11 and f DB 2 %(6+ [E+ [2) ——
C—IE‘IQ) is the dielectric boundary force density. Here, E33 = —V¢i are to be
evaluated at the interface F.

For the ionic boundary force FIB, since )- has its non-zero value only in the solvent,
A = H Jr(n — R+(I‘)). Therefore, VA contributes to a delta function evaluated at
points on the interface where the mesh lines and the interface intersect. Consequently,

the ionic boundary force is

FIB = / f’Bds (4.28)
I‘

N
= / kTZq(e—qi¢+/kT—1)ds, (4.29)
F 2'

74

where (15+ is to be evaluated at the interface P and f DB = kT%q(e—qi¢+/kT _ 1)

1
is the ionic boundary force density. At these intersecting points, the values of (6+
can be obtained by the interpolation of two grid points outside the interface and a
fictitious point inside the interface. When singular geometry occurs, interpolation
using one grid point outside the interface and two fictitious points inside can provide
an alternative. The calculation of the ionic boundary force is very similar to that of the
dielectric boundary force technically but it has less impact on the total electrostatic
solvation force. Therefore, we just mention this aspect briefly here without providing
further consideration and numerical results.

In this chapter we carry out dielectric boundary force evaluations by using both
the technique based on spline-smoothed dielectric layer [61] and our technique based
on sharp molecular surface. The former is an adaptation of the idea presented in
Refs. [61] and [45] into our code. The latter is developed in the next subsection. The
molecular surface is generated by using the MSMS [102]. A technique developed in

Ref. [128] is used to convert the triangle surface representation from the MSMS into

the required Cartesian representation.

4.2.4 The MIB solution of the Poisson Boltzmann equation

Rigorous interface technique, the MIB method, has been described in a series work
for the solution of Eq. (4.30) [42, 121, 128]. To establish notation, we present a brief
description of our MIB based PB solvers.

Assume that an interface F divides a domain Q into two disjoint parts, a molecular
subdomain $2- and a solvent sudornain 0+. The outer boundary of {2+ is chosen to be

regular. Without loss of generality, let us consider the linearized Poisson Boltzmann

 

(a) (b)
Figure 4.1: An illustration of MIB schemes
(PB) equation of the form
Nm
—V - (e(r>V4(r)) + 626246) = 44262-66 — r.)- (4.30)

1:1

where the ionic strength K.(I‘) is defined as n2 = $32 if ciqi2 The continuity con-
ditions of the solution and its flux across the interface are given by Eqs. (4.21) and
(4.22), respectively. In the MIB method. these conditions are to be implemented
rigorously in order to obtain the highly accurate and fast convergent solution to the
PB equation. Direchlet far-field boundary condition is used for Eq. (4.30).

The essence of MIBPB solvers is the treatment of interface continuity conditions
Eqs. (4.21) and (4.22) in the Cartesian grid. The standard finite difference schemes
cannot maintain the designed order of convergence due to the discontinuity in the
flux. To restore the convergence, we compute the derivatives in the PB equation
near the interface without directly using the function values across the interface.

Instead, we use only function values on the same side of the interface and their

76

smooth extensions across the interface, which are often called fictitious values, see
Fig. 4.1. A fictitious value is needed when a finite difference scheme reaches across
the interface. To determine fictitious values, we make use of continuity conditions at
the interface. This in turn enforces the continuity conditions.

The MIB method takes the dimension splitting approach which reduces a mul-
tidimensional interface problem into 1D ones. Therefore, unlike the HM, the MIB
method normally does not use multidimensional interpolation or multidimensional
Taylor expansion. Consequently, the MIB method simplifies the local topological
relation near an interface, which is crucial for 3D problems with complex interface
geometries.

Consider a case where the interface F intersects the grid in the x-direction at a
point (to, j, k), which has two nearest neighbor grid. points (2', j, k) and (2' + 1, j, k).
Fictitious values f +(i, j, k) and f ‘(i + 1, j, k) are required in a second order finite
difference scheme in the :r-direction near the intersecting point. In general, at every
intersecting point of the interface and mesh lines, there are two original interface con-
ditions and two additional first order interface conditions. These four conditions can
in principle determine four fictitious values. However, these four interface conditions
involve six first order derivatives (75:, 6);, d); , (by— , (15:, and (15;. The evaluation of
these derivatives can be very difficult for a complex geometry. In the MIB method,
we simultaneously determine as fewer fictitious values around an intersecting point
as possible so that we have the maximal flexibility in avoiding the determination of
many first order derivatives, which are often difficult to evaluate due to the geometric
constraint. Nevertheless, we have to determine at least two fictitious values so that
both of the original two interface conditions can be implemented directly or indirectly.
Consequently, we determine only two fictitious values around an intersecting point

77

and thus, use two interface conditions to eliminate two first order derivatives. The
remaining four derivatives are to be approximated using appropriate grid function
values near the intersecting point.

In order to derive two additional conditions, we define a local coordinate (E, 77, C)
such that 5 is along the normal direction, and n is in the :1: —y plane. Two coordinates

are connected by the transformation

 

 

 

 

 

 

 

F ' " "l " ' F' '
5 sin 11) cos 6 sin 11'; sin 6 cos 11) a: a:
17 = — sin 6 cos 6 0 y = P ' y - (4-31)
C — cos 1,11 cos 6 — cos 11' sin 6 sin 11') z ] z

- .- _ .1 h- . .J

 

where 6 and 1,6 are the azimuth and zenith angles with respect to the normal direction
g, respectively. Two additional continuity conditions are generated by differentiating

Eq. (4.21) along two tangential directions, 77 and C

[4,] = (412221 + 431122 + 421223) — (4.22221 + 4.71222 + 4.72223) (4.32)

[4%] = ($271931 + <15; P32 + 05:1933) - ($51931 + (15371932 + ¢§P33)- (433)

where pij is the ijth component of the transformation matrix p. From Eq. (4.31)

we have [(25, (157,, (de = p - [(1533, 65y, (DAT, where T denotes the transpose. Three flux

78

continuity conditions, Eqs. (4.22), (4.32) and (4.33) can be cast in the form

 

 

 

 

4:
’ ‘ ‘ c5;
[WE]
$3
[Q50] 2 C - , (4.34)
4’17
[$4] +
_ - ¢z
91).?
where
r . i + — + — + — .
Cl P111jl —P11/5 P123 -P12..3 P136 -P135
02 C2 = P21 -P21 P22 -P22 P23 -P23 ’ (4'35)
C3 P31 —P31 P32 -P32 P33 -P33

 

 

 

 

where C,- represents the 2th row of matrix C. As discussed above, there are six
derivatives, (bf, 66;, 6b; , at; , 6):, and 65;, in restricted subdomains. Some of these
conditions are very difficult to evaluate for complex biomolecular interface geometries.
Therefore, in the MIB method, we eliminate two derivatives as the follows. In Fig.
4.1(b) we illustrate a case where the (by. is difficult to compute and is to be eliminated.

In general, after the elimination of the lth and mth elements of the array [66}, (b; ,

79

(1'); , ¢y" , (15:, 65;], Eq. (4.34) becomes

a. [3455] + b [6,] + c [6C] = (aC1 + ng + cC3) . , (4.36)

 

 

where

a = C2zC3m - 03102771 (437)
b = C3ZC1m - 01103"-

C : CIICZm — C2lclm-

Finally, Eqs. (4.21) and (4.36) are used to determine two fictitious values near the
interface, together with the evaluation of four partial derivatives. This procedure is
systematically repeated to determine fictitious values along other mesh lines and at
other intersecting points. This dimension splitting approach effectively reduces a 3D
interface problem into 1D—like ones locally.

For the situation depicted in Fig. 4.1(b), after elimination of q"); , we compute
6);] . However, on the line passing through (3'0, j, k) in the y-direction, there is no grid
point. Therefore, 653'] has to be evaluated by using interpolation schemes with solution
values from the neighboring points. Specifically, (6,] at (4'0, j, k) will be interpolated

by using values at points (to, j — 1, k), (2'0, j, k) and (2'0, j + 1, 1;). Because these three

80

points are not normal grid points, we have to further obtain the values at (20, j — 1, k)
by interpolation using grid values at (2+ 1,j— 1, k), (2+2,j — 1, k) and (2+3,j— 1, k),
and the value at (20,j + 1, k) by using grid values at (2 + 1,j + 1, 16), (2+ 2,j + 1,16)
and (2 + 3,j + 1, It). For the value at (20,j, k), the fictitious value at (2,j, k) and grid
values at (2 + 1, j, k) and (2' + 2, j, k) will be employed for the interpolation.

There are some simple principles for the selection of two derivatives to be elimi—
nated and for the selection of the interpolation strategy. One is the feasibility of the
scheme and the other is matrix Optimization. We eliminate those derivatives that are
difficult to evaluate. For a given local geometry, we will choose to compute one of
two partial derivatives, (62+ and 6;, such that the evaluation in the .2: — 2 plane can
be easily carried out. We select an appropriate interpolation strategy such that the
resulting MIB matrix is optimally symmetric and diagonal for arbitrarily complex
solvent-molecule interfaces with geometric singularities. A detailed discussion of the
matrix optimization and 3D MIB methods for geometric singularities can be found in
Ref. [122]. A rigorous validation of MIB based PB solver, MIBPB-II, for molecular
surface singularities can be found in Ref. [121]. For a detailed treatment of charge

singularities, the reader is referred to Ref. [42].

4.2.5 Dielectric boundary force calculation

There are two features in the present work that differ from any of the previous work
in electrostatic force calculation. First, our force calculation is based on the rigorous
treatment of the PB equation with discontinuous solvent-solute interfaces. Therefore,
our evaluation of the dielectric boundary force density is to be based on Eq. (4.26),

instead of Eq. (4.18), for sharp interfaces. Second, due to the use of molecular

81

surfaces, instead of van der Waals surfaces, the distribution of reentrant surface effects
in the dielectric boundary force to individual atoms is a new issue. In the present
calculation of the dielectric boundary force, three issues are to be resolved. First,
one needs to evaluate El: in Eq. (4.26), which involves partial derivatives at the
interface. Second, one needs to numerically carry out the surface integral Eq. (4.26)
along the interface. Finally, one needs to distribute the calculated boundary forces
from reentrant surfaces to individual atoms. For the first issue, the calculation of the
surface density function is actually at the intersecting points of the interface and the
meshlines. Some of these intersecting point values are actually computed in solving
PB equation, and by modifying part of our MIB schemes, derivatives of <25 with respect
to :r,y and z~directions on these points in both inside and outside of the interface
can be conveniently computed. For the surface integration, we employ the method
developed by Smereka [106] with a minor modification. In this method, the numerical
surface integral is calculated by volume integral with integrand supported by discrete
delta functions defined only on grid points close to the interface. For the issue of
force distribution, since our numerical integral is in the form of the summation of
the product of the force density and curved small area, we can distribute the force
to individual atoms based on the location of the force density. The details of surface
density evaluation, surface integration and force distribution are given in the following

three subsections.

Evaluation of dielectric boundary force density on the interface

The dielectric boundary force on the interface as discussed in. the previous section is

given as Eq. (4.26). To compute E+ and E‘, we need to evaluate the gradient of (b

‘6

on the interface from both inside the biomolecule, denoted as —”, and outside the

82

biomolecule but in the solvent, denoted as "+”, i.e., (25.; , (bi, tag, of , ,2; and cf); .
Having obtained the solution (1) to the PB equation in the entire domain, suppose
now one wants to evaluate these partial derivatives at ($0, yj, 2k) from Fig. 4.1(b),
(b; can be interpolated by values of (b at (2 — 1, j, k) and (2', j, k), and the fictitious
value of d) at (2+ 1, j, k), which is the linear combination of values at some grid points,
while if); can be interpolated by the values of d) at (2 + 1, j, k) and (2, j, k), and the
fictitious value of at (2 — 1, j, k). In our MIB scheme, two among a5; , (15;, ab: and
d); can also be interpolated by values at one fictitious point, two grid points and six
auxiliary points (actually just regular grid points). Here, one only needs one from
these two, i.e., o; . Then, one can use Eq. (4.34) to compute 9b; , (f): and ¢;. To

this end, one rewrites Eq. (4.35) into the following form by block matrices.

 

[“196] 23;? 4527
J = [(157)] = CA - (25; + CB - 2' (4.38)
Ll¢cld _¢il_ _¢;_

 

 

 

 

 

where C A and C3 are the first three and last three columns of C, respectively. With

this setting, one can obtain the three unknown partial derivatives by

 

 

-¢,7- ( '23)
,9; =Cg1. J—CA- 4,; . (4.39)
-4537. K W-)

 

 

 

 

Note that under some rare situations, matrix CB could be singular, and one
needs to first compute four of six partial derivatives of g!) by interpolation with grid

point and fictitious values, and then use the interface jump conditions to solve the

83

remaining two. Moreover, if an intersecting point is exactly a grid point, some special
treatments are required as described in the MIB schemes [122, 123]. These treatments
are implemented in our software package and are omitted here since their description

is lengthy.

Numerical surface integration in Cartesian grids

According to Smereka, the surface integral of a density function in Cartesian grids

can be approximated by [106]

ff f(rc, y, ads = [9 fm, z>6<d(:c.y. mm z E f(:v.-.yj, zk>5.,j,kh3 (4.40)

mkk

where ($,,yj, 2k) is the coordinate of grid point (2', j, k), d(x,y, z) is the distance of
a point (:r,y,z) defined in Q from I‘, and f (any, 2) is the surface density function
defined in F. With (I) as the level set function such that (I) = 0 defines the interface,

the first order discrete delta function at (2, j, k) is given by [106]

‘ _ ~(+1?) ~(~17) ~(+31) “(11) ”(+2) “(-2)
61,j,k — dijak +6l,j.k +6l,j,k +6l.j,k +63,jyk +6i,j,kl (4.41)

84

where

 

 

 

 

 

 

cr- - D0<I>- . .
~(+:r) 2.2615”? xvez’ci’il. if (1)2,j,kq)2'+1,j,k S 0;
a,” = l x We” 0 ml (4.42)
0 otherwise
lq’i—l j kDg-‘I’z'j kl .
~ ’12 D-Q, , C, , If (Di,j,kq)i-I,j,k < 0’
diff] = | 2: i.j.k||Vo‘1’z'.j,/c|
13],.
0 otherwise
|¢ij+1k03¢ij kl
~(+ ) h2 Dicp 3 V6 L1,, if q’z‘.j.kq’i.j+1,k S 0;
0 otherwise
I‘Dz' j—1.kD39(pij,kl .
~(_ » h2 Diq, ' V41, . 1fq’i.j,kq’i,j—1,k < 0;
61.3.]; = l y 23131:” 0 23231:!
0 otherwise
I‘I’z‘ijDg‘I’ijA-l
~(+ ) h2 Di»; Ve’q,’ 1hpt,-yogi“, .<. 0;
62.3.7; 2 | z 2,j,kll 0 z',j,kl
0 otherwise
[(1). . DOq). . [
2,],k—1 z 2,],k .
" ) h2 D_ C 1 If Qiajak¢lsjak‘l < 0;
if: I z¢i,j,kllvo¢’2,j,kl
2,39
0 otherwise

where D;<I>i,j,k, US$203}: and Dg<I>i,J-,k are. upwind, downwind and central finite

difference schemes in the a-direction (a = .r,y,z), respectively, and [V5@i,j,k| =

 

V/(ng)i’j’k)2 + (D3¢i.j,k)2 + (ng)2,j.k)2 + C With C 2 1.0—10.

In our case, the function values are available only on the interface. Therefore, a

modification of the above formalism is required. After some derivation, we have the

85

following integral expression

lnxl

(f($01ijzk) h

 

+ fixayo. 20%;] + f($z‘,yj, Zoll—hE-I ’13,

(4.43)

Af(x,y,z)dSz 2

(23131061

where (rmyj, 2k) is the intersecting point of the interface and the a: meshline that
passes through (2, j, k), and n1 is the :1: component of the unit normal vector at
(11:0, yj, zk). Similar relations exist between (132-, yo, 2k) and 72y, and between (1232:, yj, 20)
and 72,3, respectively. Since Eq. (4.43) has already taken into account for the contri-
bution from irregular grid points outside the interface, the summation is restricted to
I , the set of irregular grid points inside or on the interface. Based on Eq. (4.40), the

derivation of Eq. (4.43) is quite straightforward, and is illustrated with 53:).

09:81 2' k . a1 .
’ ’ are numeric a roxrma—
V 3:.ch pp

(1) In the expression of (TS? , . .

2.1.1:

 

 

¢i+1 ',k l
ESL] and

 

tions of h; : (13,-+1 —:z:0) and [221,], which can be analytically evaluated since one
knows the interface location and normal information of the molecular surface

at its intercepting point.

(2) For grid point (2, j, k) inside the interface with a negative (I) value, if its neighbor-
ing point, e.g., (2 + 1, j, k), has a positive (I) value, the discrete delta function at
(2, j, It) will have a component corresponding to 53;) in Eq. (4.42). Meanwhile,
the discrete delta function at (2+ 1, j, k) will have a component corresponding to
52:31: in Eq. (4.42). The sum of the two components at (2, j, k) and (2+1, j, k)
, 72
IS exactly I—hfl.

(3) For a sharp interface, it is possible that both sides of (2, j, k,) have positive CI)
values, then Eq. (4.43) will have two terms in the x-direction at (xj,yj,zk)
and (IO—,yj,zk).

86

 

(a) Molecular surfaces (b) Force distribution

Figure 4.2: Dielectric boundary force calculation.

(4) In Eq. (4.40), the density function f is evaluated at grid point (xi,yj,zk). For
our application, the density function is defined only on the interface, therefore
we use interface values f(:ro,yj,zk), f(1:,-,yo,zk) and f(x,,yj,zo) instead. In
fact, the results of surface integration evaluated at the interface are better than

those evaluated at grid points, according to our numerical tests.

(5) Eq. (4.43) is also valid for the case where the interface exactly crosses a grid

point.

The distribution of dielectric boundary forces to individual atoms

In the previous two subsections, we discuss the evaluations of numerical surface inte-
gral and dielectric boundary force. To compute dielectric forces on individual atoms,
we need to distribute each unit dielectric boundary force, i.e., the force density mul—
tiplied by the associated area element, to related neighboring atoms and then collect

all the forces on each specific atom. This procedure requires to classify the interfaces

87

points into several categories since different treatments will be applied.

Fig. 4.2-(a) shows an example of the molecular surface for a small system com-
posed of 5 atoms. The red dots denote the centers of the probe, the “+” denote the
centers of atoms, and the capital letters denote the contacting points between the
atoms and the probe. The molecular surface consists of three type of facets. The
contact surfaces are parts of atomic surfaces that directly contact the probe. The
toric surfaces are the collection of the curves on the probe connecting two contact
points, when the probe rolls over two contacted atoms, as shown by the curved face
CE DB on the graph. The reentrant. surfaces are the curved facet on the probe when
it contacts more than three atoms, shown as ABC or DEF on the graph.

We first classify interface points as two types, on-atom points and off-atom points.
For molecular surfaces, an on-atom point is a point that is on a contact surface while
an off-atom point is a point that is on a toric surface and /or a reentrant surface. For
the first type, we just need to find points which are in the distance of atom radius
from the center of an atom, and attributing it to the atom. For the second type, we
need to figure out which atoms the probe contacts with. Our method is to locate
the probe center by a distance of the probe radius from the interface point in the
surface normal direction. If the distance between an atoms center and the probe
center is approximately the sum of atomic radius and probe radius, then the probe is
regarded as contacting this atom. For numerically generated MSS. a small tolerance
(threshold) is used to classify interface points. In case of the MSMS surface, the
results of the classification can be checked by those reported by the MSMS program.
Approximately, the optimal threshold is 0.013 when the MSMS density is set to 10,
and 0.0013 for density 100.

The integral method described in previous section can give a force element on a

88

small part of the molecular surface. A force element is directly attributed to an atom
if the corresponding interface point is an on-atom point. However, if an interface
point is on a toric or reentrant surface, we need to attribute the force element to all
the contacted atoms.

For a toric surface, the normal direction of a surface element is in the same plane
with directions pointed from the centers of the two atoms to the center of the probe,

the distribution of the force element to two atoms can be calculated by

F'fl—F°f2(f1-f2)

 

 

F 4.44
1 1 —(f1-f2)2 ( )
F'fz-F'f1(f2'f1)
F = f 4.4

where F is the force from the surface point to the probe center, which will be dis-
tributed into two forces F1 and F2 pointing from the two atom centers to the probe
center, respectively. Here, f, f1 and f2 are corresponding unit force vectors.

For reentrant surfaces, the probe might contact more than three atoms at a time.
These situations are dealt with as the following. First, if the probe contacts only three
atoms, the force vector is decomposited into three vectors pointing from each atomic
center to the probe center. The magnitudes of three components are computed by

solving the following linear system

f1(1) f2(1) f3i1) F1

f1(2) f2(2) f3(2) F2 =FT, (4.46)

f1(3) f2(3) f3(3) F3

- d i- d

 

 

 

 

where fi(j), 2,j = 1, 2, 3, is the jth component of the unit force vector fi, F1, F2, and

F3 are the magnitude of the distributed force.

89

If the probe contacts more than three atoms, we. need to distribute the force into
more than three directions. However, this distribution is not unique, except that all
atoms have the same distance from the probe. In the latter case, the force can be
equally distributed. Otherwise, the force will be distributed to three atoms that are
closest to the force element. Fig. 4.2-(b) illustrates this scheme. Here, the interface
point locates at point 0, which belongs to the probe surface ABCD, i.e., the reentrant
surface, when the probe contacts four atoms at points A,B, C and D. To avoid the
difficulties of decompositing force vector PO to four directions, we choose the closest
three atoms at points B, C and D, and distribute the vector PO into components
PB, PC and PD. Our numerical experiments indicate this scheme is a very good

approximation.

4.2.6 Special issues about MIBPB-III

The advantage of the MIBPB-III over the MIBPB-II is at the treatment of the singular
charges. In the MIBPB-III, instead of distributing a singular charge to its eight
neighboring grid points by the trilinear approach, the Greens function formulation
is introduced to separate the charge singularity. The electrostatic potential (25, as the
solution to the PB equation with the continuities of electrostatic potential and flux
density, is decomposed into three parts, i.e., (25 = d)* + 450 + <5. Here d)* (r) is the
Green’s function .
Nm
¢*(r) = Z —3’—— (4.47)

c’Ir—rz-l

90

solving the Poisson equation with singular charges, while ooh) is a harmonic function
in Q" satisfying the Laplace equation

V2960“) = 0 in Q“;
(4.48)

(1900') = —¢*(r) on F.

For the linearized version of PB equation, Eq. (4.7), the correction potential 43(r) =

9")(r) — fir) can be rewritten as
—V- (e(r)v_d3<r)) + 420040) = —42(r)4(r). (4.49)
subject to jump conditions
[5511‘ = 0 and lécfnlr = -l€¢3nlr = E’VW‘ + 950) - nlr- (4-50)

The MIBPB-III ensures second order convergence. A major advantage of this scheme
is that it is very accurate even if the grid size is as coarse as 1.0-1.2 A. The detail of
the MIBPB-III related calculation can be found in our earlier paper [42].

To calculate the dielectric boundary force with the MIBPB-III approach, we need
(35;, a; , 46:, (p; , a); , and (25;. Note that (if): ()7); , and 03: are equivalent to (bf, ¢+,
and 4);. However, (1);. , cf); and (I); are to be obtained with the combination of g5; ,cfiyfl
g5; (02, $3, 492, 9’0}, (1);, and 452. Here, e}, 45;“ , and (j): are analytically given, while (252.,
963 and 242 are interpolated from (150, the solution of Eq. (4.49). The evaluation of 43:,
(5; , (if: 43; , 1,7 , and (,5; is pursued in a manner similar to what we have described

in Section 4.2.5 for the corresponding evaluation in the MIPBP-II approach.

91

To calculate the reaction field forces, Eq. (4.17) can be modified as

N172
r” = 4WZQ2[V(¢O+¢)I.-=r,l (451.)

2:1

Nm
= ZFBF.
2:1

It is seen that the singular Green’s function cf)‘ does not directly contribute to the re-
action field force FRF . Therefore, the MIBPB-III approach for the reaction field force
calculation should be more accurate and stable than the MIBPB-II approach. This

chapter reports the electrostatic force evaluation based on the MIBPB-II approach.

4.3 Validation

4.3.1 Validation of the surface integration in Cartesian grids

The overall accuracy, stability and reliability of dielectric boundary force calculation
depend on the quality of PB solvers, the evaluation of six derivatives of the potential
on the interface, the numerical surface integration and the distribution of forces to
related atoms. The accuracy of MIBPB solvers has been extensively verified in the
previous papers [121] and [42]. The same MIB algorithm is modified to rigorously
evaluate the six partial derivatives at the interface. The basic numerical integration al-
gorithm in Cartesian grids was validated in [106]. However, this algorithm is modified
so that it is can be applied to our situation. A major difference is that the density
function values are defined and computed on the interface, instead of on irregular
grid points. Moreover, there is no analytical result for the surface area of molecular

surfaces in general. Therefore, we validate our Cartesian-grid surface integration al-

92

gorithm by calculating molecular surface areas of sphere, two-atom, three-atom and
eighteen-atom systems. Table 4.1 shows the validation of the surface integration with
spheres of radii 2 and 4A. Numerical results are compared with exact values. It is
seen that errors essentially depend on the grid resolution. Change in the curvature
does not affect the result very much. Approximately, second order convergence is
observed, although the method was designed as a first order one.

For two— or multi-atomic systems, a common probe radius of 1.4A is used in
molecular surface generations except specified. Our results are compared with those
of the MSMS. The locations and atomic van der Waals radii are (x y z Tvdw) =
(1.500 0.4 0.4 1.6) and (-3.000 0.0 0.0 2.0) for the two-atom system, and (2.000 0.0
0.0 2.0), (0.000 -1.0 0.0 2.0) and (0.000 1.0 0.0 2.0) for the three-atom system. For
the eighteen-atom system, these parameters are (:1: y z TvdW) = (2.0 0.9540 -0.6510
1.7), (-1.6690 0.2340 0.6650 1.7), (-0.4530 0.6870 0.4410 1.7), (0.7510 0.1480 -0.0400
1.7), (0.3930 0.8680 -1.3560 1.7), (-0.8230 1.7880 -1.1320 1.7), (-2.2840 0.2080 -1.4180
1.2), (-2.8880 1.6170 -0.4830 1.2), (-2.5270 -0.3680 0.9970 1.2), (-1.4260 0.9800 1.4350
1.2), (—0.1960 -1.1890 1.3850 1.2), (-0.7010 -1.4410 -0.3200 1.2), (1.0070 0.8940 0.7270
1.2), (1.6120 —0.5150 -0.2080 1.2), (0.1490 0.1210 -2.1260 1.2), (1.2510 1.4700 -1.6880
1.2), (-1.0810 2.2910 -2.0760 1.2), and (-0.5750 2.5430 -0.3710 1.2).

Table [4.2 gives the results of our numerical integration. It is seen from the ta-
ble that when the mesh is refined, the computed total areas converge to those of
the MSMS, which validates the numerical integration carried out at the intersecting
points of the interface and meshlines and using the corresponding normal directions.
However, the areas of contact surface, toric surface, and reentrant surface do not
Show a similar level of consistence with those of the MSMS. This is due to the fact

the MSMS surface was not analytically given but was in the triangular representa-

93

tion. We have used a threshold to categorize the MSMS molecular surface points
into contact surface, toric surface and reentrant surface. For different cases, the best
threshold may vary, which leads to errors. However, our tests indicate that the error
in classification does not substantially affect our force calculation because the pro-
jection of a surface force element into atomic center directions is insensitive to the

element classification.

4.3.2 Validation of force evaluation

Methods of validation

The validation of electrostatic force calculations is a difficult issue because of the
lack of analytical values even for small molecular systems. However, there are two
ways to evaluate electrostatic forces and to examine the consistence of force calcula-
tions. In the first approach, an atom in a molecule is moved forward (or backward)
along an axis by a small distance, e. g., 0.005A, at a time. The electrostatic free energy
of each configuration is evaluated and the electrostatic force corresponding to each
movement can be computed. For example, the force due to a small displacement (5 in
the x-direction on the 2th atom centered at (2,, y,, 21-) can be calculated as

AG(-’132:+ (5.312“, 22') - A0032 + 5. yr, 22)

F(AG) = 26 .

(4.52)

 

The other way is to obtain electrostatic solvation forces of the same movement by

summing over components

F(Sum) = ram) + F(DB) + F(IB) (4.53)

94

where F (Sum) is the total force, and F (RXN), F (DB), and F (IB) are the reaction field,
dielectric boundary, and ionic boundary forces, respectively. The consistence of the
electrostatic forces calculated in two independent approaches, F (AG) and F (Sum),
can be used to validate force schemes. As the accuracy of the MIBPB solvers has
been extensively validated [42, 121, 128], the consistence in two independent force
calculations provides a reliable check for the proposed method of electrostatic force
evaluation. Moreover, the present results can also be validated via a comparison with
those in the literature, especially those of Im, Beglov and Roux [61].

Case 1: Two-atom system

The two-atom system is a benchmark test for the validation of electrostatic force
calculations [45, 61]. Im et al. have provided a detail analysis for this system using
their PBEQ with smoothed dielectric functions [61]. In the present project, we have
reproduced their codes and results for a comparison with ours. We have set the ionic
strength to 0 so that the comparison can be carried out. The MIBPB-II is employed
in our force evaluation.

Consider two atoms of radius 2.0A. Each atom has a unit charge located at its
atomic center. One atom is fixed at (-3, 0. 0) and the other moves from (~3,0,0) to
(4,0,0). The reason for us to move the atom up to (4,0,0) is that at about (3.2,0,0),
the MS breaks up into two pieces and admits two cusps. Fig. 4.4 depicts the results of
the proposed MIBPB and of the PBEQ, which are. generated based on the algorithm
described by Im et al. [61]. The units for solvation free energy and force are kcal/mol
and kcal/mol-A, respectively. The probe radius of 1.4A is used to generated molecular
surfaces.

Detailed compasions of the results from two methods are given in Fig. 4.3. Mesh

size 12 = 0.21 is used to test these methods. The results calculated by the PBEQ with

95

 

 

     
    

kwl/mol

1%" . ‘. sag

o

kml/mol xA

 

 

 

 

 

 

 

36v
2

34 33—2—101234

 

 

kwl/mol xA
—8 N
01 O
Q
Q
O
f
'

.5

O

y
I

 

 

 

 

 

 

 

 

kcaI/mol xA

 

 

 

 

 

 

e (f

Figure 4.3: Comparison of electrostatic forces of a two-atom syitems obtained at
h = 0.21A. Red solid line: present MIB based method with molecular surfaces; Blue
dotted line: PBEQ with smooth interface. (a) Solvation free energy of the system;
(b) Total force on the moving atom calculated by using solvation free energies; (c)
Total force on the moving atom calculated from the summation of reaction force
and dielectric boundary force; ((1) Reaction force, (e) Dielectric boundary forces; (f)
Dieiectric boundary forces when the moving atom on zit-axis is in between 3.0A and
4.0

96

smoothed dielectric boundaries and by the MIBPB With sharp molecular surfaces are
illustrated with dotted blue lines and shown with solid red lines, respectively.

From Fig. 4.3(a), one sees a similar pattern of the electrostatic solvation free
energy, which confirms the reliability of two methods. We also noticed that there are
about 30 kcal/mol difference in two results, which is due to the fact that the PBEQ
uses a smoothed dielectric interface as discussed by [108, 109, 115]. The accuracy of
the MIBPB solvers has been extensively validated [42, 121, 128].

The total forces calculated by using the electrostatic free energies are given in
Figs. 4.3(b) for the PBEQ and the MIBPB. Results from two methods show different
patterns. We do not have an explanation for the oscillation of the PBEQ total force
curve.

The total forces obtained by the force density approaches are given in Figs. 4.3(e)
for the PBEQ and the MIBPB. The focre curve from the MIBPB is almost completely
smooth and monotonically decay. The smoothness indicates that the present method
is highly stable. The monotonically decay in the total force curve is consistent with
underlying physics of two separating charged spheres.

Fig. 4.4(d) gives a comparsion of the reaction field forces obtained with two meth-
ods. These forces from two methods have similar patterns. However, the peak force
value of the PBEQ is about 40% higher than that of the MIBPB. This is consistent
with what reported in the literature that the smoothed dielectric boundaries result
in 10-40% over-estimation of solvation energies and forces [108, 115].

It is difficult to obtain stable dielectric boundary forces as shown in Fig. 4.3.
Physically, we expect a smooth decay of the dielectric boundary forces as two atoms
separate from each other. The dielectric boundary forces of the MIBPB shown in
Fig. 4.3(e) have a monotonically decaying pattern, which leads to the same pattern

97

in the total force calculation. The dielectric boundary forces of the PBEQ shown in
Fig. 4.3(e) undulate over the coordinate from -3 to 1, due to the use of smoothed
dielectric functions. Moreover, the PBEQ dielectric boundary forces have oscillations
of magnitude as large as 1.5 kcal/mol. In contrast, the MIBPB results are essentially
oscillation free, even if the sharp molecular surface interface is employed.

The dielectric boundary force of the MIBPB shown in Fig. 4.3 (e) has a gradually
decreasing pattern, which leads to the same pattern in the total force calculation. Fig.
4.3 (f) gives as a zoom-in of Fig. 4.3 (e) at the end, where the two atoms are quite
away from each other. Physically, we expect the decay of the dielectric boundary
force. This pattern is observed in the results calculated by two methods. The PBEQ
result has oscillations with its magnitude as large as 1.5 kcal/mol. The MIBPB result
is relatively smoother.

It is to point out that at 1‘ == 3.2A, the MS breaks into two species, and has
two cusps. Fig. 4.3 (f) gives as a zoom-in of Fig. 4.3 (e) at the end, where the two
atoms are quite away from each other. As expected, our MIBPB forces (F (Sum)) Show
no problem for this geometric singularity at both mesh sizes. Note that undesirable
spikes are from F (AG), rather than from (F (Sum)) The PBEQ result has oscillations
with its magnitude as large as 1.5 kcal/mol.

Finally, Fig. 4.4 gives a summary of results calculated by the PBEQ and MIBPB
at two different grid sizes. The PBEQ results with smooth dielectric boundaries at
h = 0.21A, as shown in Fig. 4.4(a), are identical to those of Ref. [61], which validates
our PBEQ code. The PBEQ results at h = 0.21A indicate a good consistence of
F(AG’) and F (Sum), which validates the PBEQ approach. Similarly, MIBPB results
at h = 0.21A, as shown in Fig. 4.4(c), demonstrate a good consistence of F (AG) and

F (Sum), which indicates the reliability of proposed method. For both methods, less

98

 

    
  

 

 

 

 

 

 

---.F (A G)
- - -F (RXN)
3° - - . F (DB)
— F (Sum)
20
E
E 10
:8; :’ i.
0 ' ’ ". u.‘;.“.".’\"l\"5'¢".'.fll"
0"!” . l i,
-10.
-2 4
93 '2 -1 0 2 3 4
(a)
40 — F (A G)
- - - F (RXN)
- - F (DB)
30 — F (Sum)

kwl/mol xA

 

-1G

 

 

 

 

 

-3 -2 -1 0

n

(v)

kcal/mol xA

 

—-F(AG)

- --F(RXN)

. . I-v-'F(DB)
I ' —F(Sum)

 

 

 

 

 

l
40 —F (A G)
- - - F (RXN)
- - F (DB)
30 ' -——F(Sum)
20
10
o "y w ‘
0‘4 0 o 1“ Q.“v"\a“""'.-
-1O

 

 

 

 

 

 

 

—3-2-161234

(d)

Figure 4.4: Electrostatic forces calculated for two separating charged atoms. The
vertical axis is the position of the second atom and the horizontal is the forces. (a)
PBEQ at h = 0.21/3.; (0) PBEQ, at h = 0.42131; (b) MIBPB at h = 0.214; (d) MIBPB

at h = 0.424.

40'

 

     
  

 

 

 

 

_ F (Tot FD)
30 - - - F (RXN)
....... F (Contact)
20 .---‘ F (Toric+Reen) - “3'“:
-— F (Tot Sum) ’ ’ I
§ 10 o r p U. ’ ”
6 .1" I." "" PM» .
E 0 v"-.. I, ~_—' “‘.'.“"5'
z .Ii
.2 -10.
-20
-30 .
-40 ‘ l I I
-2 -1 0 1 2

Figure 4.5: Forces calculated for a three-atom system

numerical fluctuation is observed in the summation of force components.

It is interesting to examine the ability of the proposed MIBPB method in using a
relatively large mesh size. Figs. 4.4(b) and ((1) present the force calculations obtained
at the mesh size of 0.42A. The PBEQ predicts reasonable reaction field forces and its
F (AG) and F (Sum) show a good consistence. However, the dielectric boundary forces
of the PBEQ suffer from dramatical oscillations, indicating the approach is unstable
and unreliable at this mesh size. Note that the instability in F(AG) means that the
electrostatic free energies obtained with the PBEQ are unstable. The MIBPB-II is
used in this test too. Its F (Sum) as shown in green color is still very smooth, and
shows a good consistence with that computed at the mesh size of 0.21A.

Case 2: Three-atom system]

Since molecular surfaces admit reentrant surfaces, it is important to further vali-
date our algorithm for a 3-atom system. In this case, one atom is moved from (-2,0,0)
to (2,0,0) in the 113-direction while the other two atoms are fixed at (0,-1,0) and (0,1,0).

All atoms have radius 2.0A and charge 1.0ec. From Fig. 4.5, it is seen a good consis-

100

tence of the :r-direction forces on the moving atom obtained by two different methods.
The symmetric pattern about :r. = 0 confirms the symmetric setup of the testing sys-
tem. One thing that is worth of mentioning here is that the forces calculated by
electrostatic free energies show more oscillations, which in fact depend on the quality
of the MSMS surface generation. These oscillations can be avoided or created with
different setting of the MSMS parameters. In contrast, the total forces calculated
by the summation of force components are much smoother and more stable. This is
important for PB based molecular dynamics simulations.
Case 3: Eighteen-atom system

The atomic coordinates (any. :5) and van der waals radius (1‘) for the eighteen-
atom molecule Cyclohexane are given as (:r, y, z, r): (:51 0.9540 -0.6510 1.7), (-1.6690
0.2340 0.6650 1.7), (-0.4530 -0.6870 0.4410 1.7), (0.7510 0.1480 -0.0400 1.7), (0.3930
0.8680 -1.3560 1.7), (-0.8230 1.7880 -1.1320 1.7), (-2.2840 0.2080 -1.4180 1.2), (-2.8880
1.6170 -0.4830 1.2), (~2.5270 -0.3680 0.9970 1.2), (-1.4260 0.9800 1.4350 1.2), (-0.1960
-1.1890 1.3850 1.2), (0.7010 -1.4410 -0.3200 1.2). (1.0070 0.8940 0.7270 1.2), (1.6120
-0.5150 —0.2080 1.2), (0.1490 0.1210 -2.1260 1.2), (1.2510 1.4700 -1.68801.2), (-1.0810
2.2910 -2.0760 1.2), and (-0.5750 2.5430 -0.3710 1.2). The partial charges of these
atoms are artificially assigned to the atomic centers as 1 cc. In this test, we fixed the
coordinates of atom 2-18 while move the :r-coordinate of the first atom 2:1 from 1 to
3. In Fig. 46(b), electrostatic forces on the first atom obtained by the two different
methods show a good consistence.
Case 4: A protein molecule (1ajj) with 519 atoms

The final example is chosen as the molecular surface of Homo sapiens (PDB:
1ajj) with 37 residues and 519 atoms. Fig. 4.7(a) shows its structure. We move
the outmost nitrogen natively located at (—0.169,7.698,13.415) along the positive and

101

 

 

 

 

 

 

 

 

 

—--F (A G)
400 ' ' 'F (RXN)
F (Contact)
- - F (Toric + Reen)
300 . —F (Sum)
200» _ . — - ' '
§ ' v ’
6 ..:..,. a ’ ‘ ‘
ti ‘°° Mv
.2 _¢“ ~l~lo:lh "'2
0 * ~ “~., ‘
"""1 1.5 2 2.5 3
(a) Molecular Surfaces (b) Forces

Figure 4.6: An eighteen-atom molecule Cyclohexane

2.
.3...-
I.
~'\
0 ~.--
§
--

  
  
  

I‘ll]:--------------:>I--..‘.___

 

 

 

 

 

 

_2 . ‘
E -4- """" . ..
0 .fi,
5 2
(U u
e —- F (A G)
-3 - - - F (RXN)
------ F (Contact)
—10 - - F(Toric+ Reen)
—— F (Sum)
-12 A + ._ . m l l . A
-O.2 -0.1 0.0 0.1 0.2
(a) Molecular Surfaces (b) Forces

Figure 4.7: Homo sapiens

negative direction of the .r-axis for a distance of 0.25A. This movement cannot be too
large, or it will violate the physical limitation. The calculation results are shown in
Fig. 4.7(a) where the confirmation of results from two methods has been achieved.
Finally, we provide an alternative way to confirm the convergence of the resulting
forces from our method. In Fig. 4.8, we plot the correlation of forces obtained at
different mesh sizes, i.e., 0.25A, 0.5A and 1.0A. All the forces in r, y and 2 directions

for each atom are plot here. Born the figure, we can see a good convergence of reaction

102

field forces at mesh sizes of 1.0A and 0.5A compared with that obtained at 0.25A.
we expect that the MIBPB-III based force correlation will have a better performance.
For dielectric boundary forces, as well as their toric and reentrant components, their
accuracy is significantly reduced at the mesh size of 1.0A, because of the error in the
surface integration at such a coarse mesh. However, at 0.5A, the convergence is good
as seen in Figs. 4.8 (c)-(e).

Similar correlation diagrams are given for PBEQ forces as plotted in Fig. 4.9.
Although reaction field forces obtained at the mesh sizes of 0.25 and 0.5A are shown a
good consistence, their correlations between 0.25 and 1.0A are quite poor. Moreover,
the dielectric boundary forces calculated at the mesh sizes of 0.25 and 0.5A differ
dramatically. Based on these results and the results shown Fig. 4.4, we call for the
attention that there is a further need to reexamine the convergence, stability and
reliability of previous PB based MD methods that have not done the convergence

studied.

4.4 Application to molecular dynamics simulation

4.4.1 The implementation of MIBPB based PB forces in the
AMBER package for molecular dynamics

The present MIBPB calculates the electrostatic potential and provides electrostatic
forces asserted on each individual atoms for solvated biomolecules. However, to per-
form the molecular dynamics simulation, it takes more computational modules such
as other components of the force field, initialization, randomness, time evolution and

so on. In this work, we implement the MIBPB 1n the AMBER package. Specif-

103

Toric Force at 0.5A RXN Force at 0.5 A

DB Force at 0.5A

 

30

20

10

-20*

 

 

 

“3930 ~20 ~10 0 10 2‘0 30

15

RXN Force at 0.25 A

(a)

 

10

-10

 

 

A

 

-1 x r A g,
-§l5 -10 -5 0 5 1O 15

30

Toric Force at 0.25A

(C)

 

20

10

-10+

 

 

 

 

-3 . - - . -
£50 -20 -10 0 10 20 30

DB Force at 0.25A

(e)

RXN Force at 1 .0A

Reenh’ant Force at 0.5A

30

20'

 

 

O

 

 

 

”2%

30

0 ~20 ~10 0 10 2o 30

RXN Force at 0.25 A

(b)

 

20'

10*

 

 

 

1‘20 ~20 ~10 0

10 20 30
Reentrant Force at 0.25A

((0

Figure 4.8: Correlations of forces obtained at different mesh sizes by the MIBPB
for protein 1ajj. (a) Reaction field forces between mesh sizes of 0.25 and 0.5A; (b)
Reaction field forces between mesh sizes of 0.25 and 1.0A; (c) Toric forces between
mesh sizes of 0.25 and 0.5A; ((1) reentrant forces between mesh sizes of 0.25 and 0.5A;
(e) Dielectric boundary forces between mesh sizes of 0.25 and 0.5A.

104

 

 

   

 

 

 

 

 

 

 

O
< 30 < 30*
O Q
‘0. 20 1 C3. 20* *
O I: (o O
.5 10’ i (a 10 (2 ]
g 0 1 g 0.
o u.
u.
.27 -10 E -10 1
.2 ~20 2 ~20» , l
Q 30 '2 , ” E 30 (I
g 4 £2
'40 z; '40’
")0
-5%“‘+ *‘ r -5%%W
- 0-40-30-20-10 0 10 20 30 40 - 0-40-30-20-10 0 1O 20 30 40
Recational Force at 0.25A Reactional Force at 0.25A
(a) (b)
50 . ~
< o
‘0. o _
g 25 o 0" U
8 f0 Goff?”
9 0 1.-
LL ‘3}
E Cw A
g -25 \ 4 ‘
1: .
8
m -50
U Q
'1:
E) '75 1
.‘2 o
D

 

 

 

-10 A g A .
-900 -75 -50 -25 0 25 50
Dielectric Boundary Force at 0.25A

(C)
Figure 4.9: Correlations of forces obtained at different mesh sizes by the PBEQ
for protein 1ajj. (a) Reaction field forces between mesh sizes of 0.25 and 0.5A; (b)
Reaction field forces between mesh sizes of 0.25 and 1.0A; (c) Dielectric boundary
forces between mesh sizes of 0.25 and 0.5A.

105

ically, we replace the original PB solver and force driver in the AMBER [85, 86]
with our MIBPB solver and MIBPB force driver, respectively. We use 2gb = 8, a
value currently not being used by AMBER-GB modules, to switch the implicit sol-
vent model to the MIBPB. Seven files in the AMBER molecule dynamics module of
Sander, such as dynlib.f, egb.f, force.f, mdread.f, printef, runmd.f and sander.f are,
to some extend, modified to meet our needs and the original functions of the AM-
BER are kept intact. We add seventeen MIBPB files, prefixed with “mibpb”, i.e.,
mibpb_bicggt.f, mibpb_fish.f, mibpb.matrix.f, mibpbsharppatchf, mibpb_chgdist.f,
mibpb.force.f, mibpb.module.f, mibpb_topology.f, mibpb.f, mibprnsmsmibw13.f, mibpru.f,
mibpb_para.f, mibpb.fftpack.f, mibpb.main.f, mibpb_poisson_solver.f, mibpb.fdweights.f,
and mibpb_kirkwood.f. Meanwhile, we modified the Makefile of the Sander to include
the compile of all the newly—added files. For AMBER users, to add the MIBPB mod-
ule, all they need to do is to replace the above mentioned 7 files and the Makefile
with our files, and add 17 new files under “\$AMBERHOME\src\sander\”.

In the mibpb.f, the file which connects the AMBER and the MIBPB, the coordi~
nates, charges, and radius of all atoms, together with some parameters, are obtained
as the input of the MIBPB. After calling the subroutine mibpb(), the molecular sur-
faces are generated, the PB equation is solved, and the solvation energy and forces
are calculated. Eventually, energy and forces are assigned to the corresponding global
variables in the AMBER, e.g. (eelrf for solvation energy and pbfrc for forces), as the
output of the MIBPB. The cut-in place of the MIBPB, i.e, the calling of the sub-
routine named mibpb(), is included in the subroutine pb_force() in the file pb.force.f,
paralleled to the calling of the subroutine pb_fdfrc.

It is also necessary for us to point out that the MIBPB, as still in the early stage of

its development, does not provide many options for the user. The molecular surface

106

is obtained by running the MSMS [102]. Charges are treated by trilinear interpola—
tion. The F F T is used to solve the PB equation in vacuum, and a preconditioned
bicongigate method is used to solve PB equation in solution.

In summary, to run the MIBPB module embedded in the AMBER, we change the
algorithm in providing reaction filed forces, dielectric boundary forces, ionic boundary
forces, and electrostatic solvation energy while keep the force fields of bond, angle,

dihedral, van der Waals, and Coulomb intact.

4.4.2 Molecular dynamics simulations

In this subsection, we demonstrate MD simulations using the MIBPB based PB ap-
proach in conjugation with the AMBER package. For a comparison, the original
AMBER PB module developed by Luo and his coworkers [85, 86] is used for the same
simulations. At mesh size of 0.3A, two small biomolecular systems are considered in
our studied. The first system is the alaning dipeptide with 22 atoms. As suggested by
the AMBER tutorial, we started from a pdb file with the first atoms (C,C,N) of the
three residue/ terminal (ACE, ALA, NME), and used the leap to add other 19 atoms.
We relaxed the system with the graphic tools provided by the leap and minimized
it with 200 steps of steepest descent runs and 300 steps of conjugate gradient runs.
After that, we run 10ps MD (10000 steps) simulations. The results are given as in
Fig. 4.10. The RMSD of the MIBPB is slightly lower that that of the AMBER PB
as shown in Fig. 4.10 (a). In Fig. 4.10 (b), it is seen that the solvation free energies
of the MIBPB are more stable and higher than those of the AMBER PB. The total
potential energies depicted in Fig. 4.10 (c) show a tendency similar to those found in

the solvation free energies. For temperature shown in Fig. 4.10 (d), since the system

107

RMSD EPB

 

 

 

 

 

 

 

 

0.5 ,3
.1» .42.).1111 21 21’ 1,; -~
03 , ‘ l[ll)llll[l [WI lll[1.1, A}, Ill (ll A
'11,)“ , |l] I] I 3 ll‘lrlfl' All i ,1»
A ‘5 ' l]1lllll ,[Anlm' Ml ], 5—2o]l l], W N?“ "W “vi", \l l/ ll M“
0.2 l ll [ll ill illlll W 32 'J 1
l] ‘21
01 ‘25
' «- MIBPB
n -—AMBERPB
”o 2 4 ps 6 a 10
(a)
EPTOT

 

 

keel/mot

    

 

 

 

'40 ] ---~ -— MIBPB ~~~~~~ ,_ MIBPB
—~— AMBER PB —— AMBER PB
45 c ,
o 2 4 ps 6 8 10 0 2 4 ps 6 8 10

 

 

 

(C) (C1
Figure 4.10: A comparison of 10 ps MD simulations of the 3.1811ng dipeptide ob-
tained with PB based methods. Results of the MIBPB and AMBER PB are illustrated
in red and blue lines, respectively. (a) RMSD; (b) Solvation free energies; (c) Total
potential energies; (d) Temperature.

is very small, the result has a limited significance.

The other biomolecular system used is the WW domain binding protein-1, the
peptide bound to 65 KDA yes-associated protein. The structure is obtained from
the complex (ljmq) and has 10 residues (GTPPPPYTVG), 138 atoms. We follow
the AMBER procedure to obtain the pdb file and add hydrogens by using leap. The
minimization and MD simulation have been run in the same setup as that of the first
system. The translational and rotational motions are removed after every 100 steps

of simulation. From Fig. 4.11, we can see that for both methods, RMSDs shown in

108

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.5 . . . . 0 1 . r
‘ ' II“; ' .3 If l {Av-4"?!“ 1 HI (I ,H ] f o .9 A”
A 1 1 1 1'11 .1 s -1511
0.2 . ' 1 “#51191 g _. "—1MIBPB
-2001 '—AMBER PB 3
h31AM. fix. ’1 A 1 A" ”"" 1‘ ft:-
0.1 . .. .. .. .. MIBPB _250 . \/,/\,\ j 'P‘Lj 1,- Ver V’ w i V\ F ,1
33 \, , WW
= ~~ —~ AMBER PB
0 ‘ ‘ 1 1 -3oo 1 1 1 1
o 2 4 ps 6 8 1o 0 2 4 ps 6 8 10
(a) (b)
EPTOT TEMP
-50 . . - . 500 1 . .
, 400 ' I ‘ 1
-100 F .1 I}. "11 '5',ng I , f, ‘ 3.1,;th , .‘Zi H WA“ ‘1’ (,fifkul‘l
("RIP 1": .f' ;'r.1 / ; if; ‘f- ., 1.,- if N -. ,- 3.11:]? lI \‘ié . ,1 g . Q , (fl ‘ :1 A.
3 ll l1” ‘ "131%., 1‘ j, 11 I"I 1 111., ‘ l .11"... 1‘ 3001 ll» ‘11 ‘1“ l” 11 .fifii"1l‘lt‘1"1f“i1€1.t‘1~’l"1i ‘-,."“1‘
I; 'J ’ ‘ U" 1 0,11 I , 1 '-
g '150 flu?” WV U ' I‘fl‘Lfif, ‘ K i. ‘
3‘2 .1 ‘1... 200 1
I '”‘ 1
-2001 ‘ 100
i . .. , MIBPB -—— MIBPB
AMBER PB 0 , . f— AMQER PB
'250 m ‘ ‘ r o 2 4 6 8 10
o 2 4 ps 6 8 10 W

(C) (d
Figure 4.11: A comparison of 10 ps MD simulations of the WW) domain binding
protein-1 obtained with PB based methods. Results of the MIBPB and AMBER
PB are illustrated in red and blue lines, respectively. (a) RMSD; (b) Solvation free
energy; (0) Total potential energy; (d) Temperature

Fig. 4.11 (a) are quite low and temperature shown in Fig. 4.11 (d) is very stable with
respect to time. There is a significant difference in the solvation free energies obtained
by the two methods. The results given by the AMBER PB indicate only oscillations
while those given by MIBPB showed a pattern of solvation energy increasing in the
first four thousand steps. Most of ten residues in the chain are non-polar, and a few
of them are polar. There is no charged residue. This kind of chains in the water will
tend to fold due to the hydrophobicity. After examining the trajectory in VMD, it is
found that the structures obtained by the MIBPB stay mostly in folded morphologies

while those of the AMBER PB exhibit mostly unfolded patterns.

109

 

4.5 Concluding remarks

Multiscale implicit solvent models have attracted much attention in the past decade
due to its efficiency in reducing computational cost by modeling the solvation and
mobile ions as dielectric continuum. The Poisson-Boltzmann (PB) equation is a
primary implicit solvent model and has been widely applied in molecular biology and
material sciences. PB based molecular dynamics (MD) methods have a potential to
tackle large biomolecules and have been extensively studied in the past ten years.
However, no interface based technique has been applied to any of the existing PB
based MD methods. The present work presents the first development of an interface
technique based MD method via the PB equation formalism. The matched interface
and boundary (MIB) method [122, 123, 126, 129, 130] is developed to handle the
sharp molecular surface (MS) in the PB based MD approach. This work is a natural
extension of our previous endeavor on the deveIOpment of MIB based PB solvers
{42,121,128}

An MIB based PB solver, MIBPB-II [121], is employed to obtain electrostatic
potentials. To calculate the electrostatic forces associated the use of sharp MSs,
new force density expressions are derived by modifying previous results of Gilson et
a1. [45]. To compute dielectric boundary forces, a complete set of potential gradi-
ents are required at all the intersecting points of the interface and all the meshlines.
These gradients are calculated by extending our earlier MIBPB scheme. The surface
integration in the Cartesian grid is a crucial issue in the dielectric boundary force
calculation that often limits the accuracy and convergence. To this end, an advanced
surface integration technique developed by Smereka [106] is adopted with a minor

modification. Unlike the van der Waals surface where each facet belongs uniquely to

110

an atom, the MS admits reentrant surfaces which do not uniquely belong to a single
atom. Therefore, forces from the reentrant surfaces have to be assigned to related
atoms. A set of force distribution schemes are introduced. The resulting electrostatic
forces due to the displacements of atoms are carefully validated by a comparison with
those obtained by differentiating the electrostatic free energy. A large test set, in-
cluding a diatomic system, a triatomic system, a cyclohexane, and a homo sapiens,
is utilized in our validation. Correlation of results obtained at different mesh sizes
are plotted to verify the convergence of the present method. The correlation patterns
indicate good convergence of reaction field forces at mesh sizes as coarse as 1.0A.
For dielectric boundary forces, good convergence can be obtained at the mesh size of
near 0.5A. The MIBPB solver and force driver are then interfaced with the AMBER
MD package, Sander, parallelizing to its original finite difference based PB equation
solver [85, 86]. The MD results obtained by the present method and the AMBER PB
method are reported for two small biomolecular systems to demonstrate the stabilities
and convergence of the present MIBPB based MD method.

A comparison of electrostatic force calculations is given to the CHARMM PB
based MD method [61]. For this purpose, the PBEQ code is regenerated in the
present work. The PBEQ [61] method relies on a surface smoothing technique to
maintain stability. The present study indicates that the proposed MIBPB based
MD method exhibits much better accuracy and stability than the PBEQ does. The
large deviations in dielectric boundary forces obtained by the PBEQ over two sets of
meshes, i.e., 0.25 and 0.5A, call for serious attention to the convergence, stability and
reliability of previous PB based MD methods that have not done similar convergence

analysis.

111

 

$0.5m momdm 38am

 

mmA mod mvom mug vod mmdm mmwod

wna 3d madam and 26 3.3 mafia

om.m mvd wmdom EA cod hodv mmd

Sufi ovd “.9me pom mo.m mmwv md
. mfiw 3&2 . om.w 2...? H

 

EEO nohm awn/x EEO Sta 83x :
< v 36mm <m .04:me :32

 

 

 

moomtsm Eoioifi now Echo :osfiwofi: "HA... 2an

112

 

 

 

mm: 8.3 5.3 $42 Sac may: 5.3. mag cod 3.2 3.3. 23% mEmE

3.3 muwv Emma 3.3: 3A 3.3 omwu 8.3 cod Sufi 3.: anew made

ombfi mfimv 8.3 8.2: EA ww.3 Each ww.ma cod mm.mH 2.2. mw.mw mad

3.3 8.3 give monH mod 3.2 5.3 mwdm cod 92$ weds mfmw md

com: 023 .38 :33 com: crop :58 :38 so“: 2.8» 4:80 :83 e
Eofifwa EopduonB 803.038

 

 

 

 

.mEmE 23 was €222ng 3288:: Bob 8333.8 398 media we semthEoO "N6. oififi

113

Chapter 5

Summary of Thesis Contribution

and Future Work

5.1 Thesis contribution

The main contributions of this thesis are the follows: the construction of a second
order convergent Poisson Boltzmann equation (PBE) solver using the Green’s function
formalism for treating geometric and charge singularies, and the development of an
interface technique based molecule dynamics (MD) method via the PBE formalism.

The efficiency and accuracy of current PBE solvers are greatly hindered by the
challenges in enforcing interface jump conditions, treating of geometric singularities
and regularizing charge singularities. In earlier work, an interface technique based
PBE solver was developed in Professor Wei’s group by using the matched interface
and boundary (MIB) method, which explicitly enforces the flux jump condition at
the solvent—solute interfaces and leads to highly accurate bimolecular electrostatics in
continuum electric environments. However, such a PBE solver, denoted as MIBPB-I

[128], cannot maintain the designed second order convergence whenever there are geo-

114

metric singularities, such as cusps and self-intersecting surfaces. Moreover, the matrix
of the MIBPB-I is not optimally symmetrical, resulting in convergence difficulty. The
second generation of MIB based PB solver, the MIBPB-II [121], was constructed
to handle geometric singularities and reduce the condition number of the MIB ma—
trix. However, when the mesh size approaches half of the van der Waals radius, the
MIBPB—II cannot maintain its designed accuracy because the grid points that carry
the interface information overlap with those that carry distributed singular charges.
The first main contribution of this thesis is to tackle this difficulty by using a Green’s
function formalism. The solution is decomposited inter regular and irregular parts so
that the solution to the singular charges is obtained in an analytical form, resulting
in extra interface flux jump conditions. The latter are treated on an equal footing
as the original interface flux jump conditions in our MIB framework. This method,
denoted as MIBPB—III [42], takes care of both geometric and charge singularities,
and is able to provide highly accurate electrostatic potentials at a mesh as coarse as
1.2A for proteins. Consequently, the MIBPB-III is a few orders of magnitude more
accurate at a given mesh size, and about three times faster at a given accuracy than
the APBS, a recent multigrid PB solver.

In their multiscale setting, PB based MD approaches have a potential to tackle
large biological systems. Obstacles that hinder the current development of PB based
MD methods are the concerns in accuracy, stability, efficiency and reliability due to
the presence of discontinuous dielectric coefficients at complex solvent-solute inter-
face. The latter has geometric singularities as represented by molecular surfaces or
many other popular bimolecular surfaces. Recently, the MIB method has been uti-
lized to develop the first second-order accuracy PB solver that is numerically stable
in dealing with discontinuous dielectric coefficients, complex geometric singularities

115

and singular source charges. As the second main contribution of this thesis, the first
interface technique based MD method is deve10ped via the PB formalism [41]. The
MIB method is used to handle sharp molecular surfaces and associated geometric sin—
gularities. The classical formulation of electrostatic forces is modified to allow the use
of sharp molecular surfaces. Accurate reaction field forces are obtained by directly
differentiating the electrostatic potential. Dielectric boundary forces are evaluated
at the solvent-solute interface using an accurate Cartesian-grid surface integration
technique. The electrostatic forces located at reentrant surfaces are appropriately
assigned to the involved atoms. Extensive numerical tests are carried out to validate
the accuracy and stability of the present electrostatic force calculation. The resulting
electrostatic forces are compared with those in the literature. The present MIB and
BP based molecular dynamics simulations of biomolecules are demonstrated in conju-
gation with the AMBER package. We show that there is a pressing need to examine

the convergence, stability and reliability of previous PB based MD methods.

5.2 Future work

The MIBPB, as a typical application of interface based methods in solving elliptic
equations with discrete coefficients and singular sources, has been developed in the
past four years in our group, connecting the PhD theses of Dr. Zhou [127], Dr. Yu
[120] and the present work. However, there are still many aspects for which developed
methods can be improved and new research topics can be introduced.

From the modeling point of view, an immediate extension is to model bio—systems
in multiple domains (number of domains greater than 2), which will require new

schemes that are different from the established method. A possible application of

116

multiple domain methods is ion-channels, which have been being under development
in our group.

All the discussions in this thesis are based on linear elliptic equations or the linear
PB equation. The extension of current approaches to cope with nonlinearity is also
needed.

As our first attempt, the present MIBPB based MD package is far from mature
and can be improved in a number of ways. First, the MIBPB-III method [42] can be
implemented which should provide more accurate electrostatic potentials at coarse
meshes. Second, the current numerical surface integral scheme is of first order in
convergence and can be improved. The accuracy of dielectric boundary forces is a
limiting fact at large grid sizes in the present method. Third, the M83 are generated
in triangular meshes by using the MSMS program [102] and converted to Cartesian
meshes by a scheme developed in [128]. This setup is very stable during the MD
simulation but does cause a stability problem on a very rare basis. There is a trade
off between the efficiency and accuracy in the choice of the surface generating density.
More optimized algorithm is expected to resolving this issue. Finally, the precondi-
tioned bi—conjugate solver is used in solving the PB equation and its matrix system
has been optimized in MIBPB-II. However, the choice of a precondition matrix and
even the decision about which linear system solver is to be used is still worthy of
further exploration. Due to these reasons, the present MIBPB based MD package
can only meet the speed requirement of small-scale systems in relatively short simu—
lation time. Various improvements of the present MIBPB based MD package are to
be carried out.

Various applications, such calculating solvation free energies, pKa values and elec-

trostatic analysis of interaction surfaces can be proceeded with the further develop-

117

ment. of the MIBPB package, especially for solvated bio-systems with intense com-
putational demand and accurate request, e.g., the chromatin folding. The MIBPB
package is expected to be online in 2008 and updates are to be provided according to

the user’s feedback.

118

Appendix

Kirkwood’s model for dielectric sphere

Kirkwood’s dielectric sphere provides an excellent benchmark for testing Poisson
Boltzmann (PB) solvers in terms of accuracy, speed of convergence and efficiency.
In his paper [67], Kirkwood provided very detailed physical analysis and solutions
to the PB equation of a spherical molecule in solvent. Here is a brief summary of
Kirkwood’s results used in the present work.

The model is a spherical molecule with dielectric boundary (radius a) and ionic
boundary (radius b) in solvent. The dielectric constant is 6‘ inside dielectric boundary
and 6+ in the exclusion layer and the solvent domain. The Laplace equation in sphere
coordinates is given by:

. 2
$53; (73%;?) + 3:13;; (smog?) + mg; = 0, (5.1)
where r is the radial coordinate, 6 E [0,7r] is the polar angle from the z-axis and

(,0 E [0, 2gb) is the azimuthal angle in the :1: — y plane with respect to the :r-axis.

The solution of Eq. (5.1) has the following expansions

00 n
(f) = Z Z (Kmnr_n_1 + Jmnr")P,T,”(cos 6)eim'*0, (5.2)

n=0 m=—n

where HT are the associated Legendre function of the first kind, and Kmn and Jmn

are constants to be determined by boundary or interface conditions.

119

0 Inside the spherical molecule, where singular charges exist, we have
Nm
—V%‘ = 47r Z q,6(r — r,). (5.3)
i=1
The potential (25‘, as the solution to Eq. (5.3), can be written as
qb‘ = (15“ + 43, (5.4)

where (15* is the Green’s function

Nm, q_
'* = —-z—— 5.5
¢ gé—lr’ril ( )

where Nm is number of charges in the spherical molecule, q,- and r, are the
charge and position of the ith atom, respectively. The Green’s function can be

further expanded as

oo Tl E .
(15* = Z Z FrinT-l—nfiPg‘(c086)ezm‘p, (5.6)
n:0m--—n
where
(n " lml)! Nm n m imcp'
Emn 2 W 291% Pn (003908 f, (5-7)
’ £21

with (r,,6,-, 90,-) being the polar coordinates of the position of the ith charge.
However, in practical calculations, this expansion converges slowly for small 1‘.

Therefore, it is more convenient to directly evaluate Eq. (5.5).

120

The non-singular term, (i), is given by

00 n
(5: Z Z (anrn)P,T,n(cos 6)eim‘p. (5.8)
n=0 m=-n
In the spherical shell, i.e. the ion exclusion layer bounded by the spherical
surfaces of radii a, and b, the potential (be governed by the equation V2¢e = 0

is given by:

00 Tl

(be = Z Z (Cmnr_"_l + Gmnrn)P,T,n(cos 6)eim‘p. (5.9)

71:20 m=—n
However, this part was not used in our model.

Outside of the sphere of radius b, the potential 05+ is governed by
v2¢+ — n2¢+ = 0, (5.10)

where n is the ionic strength defined previously, (25+ is given as

00 n

qb+ = Z Z (AMT—"‘1)e—mKn(nr)P;,n(cos0)eim‘p, (5.11)

n=0 mz—n

where
n

2371! 211 — s l s
K7103) : Z s!(2n()!(n _3)!:1: ' (5'12)

 

s=0

All the coefficients Am”, Cmn, Gmn and Bnm are to be determined by boundary

-— 6
conditions and interface condition: ct)‘ = (be, Egg,— 2 @9557, (158 = ¢+,

e +
03¢: = 59%— and ¢+(oo) == 0.

In fact, we used a two—domain model in this work — the ion exclusion layer was

121

not considered. Assuming weak ionic strength, the potential inside the molecule, gb",

with Enm defined as in Eq. (5.7), is given by:

00 Tl

45— = Z Z FELT}- + an7‘")P,T(cos 0)eim‘p, (5.13)

n=0 m=—n

and the potential outside the molecule, ¢+, is given by

00 71
5+ = Z Z (Cmnr—"-1 + G,m,r")1>,§n(cosa)eimv. (5.14)

7:20 m=—n

All coefficients an, Cmn and Gmn can be obtained via boundary and interface
conditions: (b‘ = (N, (592%: = Nag}: and ¢+(oo) = 0.

Kirkwood’s solution provides an effective benchmark for testing our scheme. As
mentioned before, the accuracy of the entire scheme for (:5 depends on the accuracy of
$0 solved from the boundary value problem, of the flux computed by differentiating
950, and of (.5 solved from the interface problem. In Kirkwood’s solution, ¢* in Eq.
(5.6) is corresponding to qb“ in Eq. (4.47). While (f) in Eq. (5.8) is corresponding to
the sum of ¢0 and <23 in (3.1). Finally, W" in Eq. (5.14) is corresponding to the (,5 in
Eq. (4.49) for the potential outside of the molecule. Meanwhile, since (1)0 satisfies Eq.

(4.48), it will have the following analytical expansion with spherical harmonics

00 7?-
$0 = Z Z (Dmn'rn)P,',"’(cos 6)ei‘m“9, (5.15)

1120 m=-n

Enm

where Dnm = W. The existence of the analytic solution of qb'“, $0, and (13 makes
6
it possible for us to test the accuracy of our scheme step by step in a sphere with

multiple charges.

122

References

[1]

[2]

[3]

[4]

l5]

l6]

[7]

[8]

[9]

[10]
[11]

[12]

[13]

L. Adams and Z.L. Li. The immersed interface / multigrid methods for interface
problems. SIAM J. Sci. Comput., 24:463—479, 2002.

H. Ben Ameur, M. Burger, and B. Hack]. Level set methods for geometric
inverse problems in linear elasticity. Inverse Problems, 20:673—696, 2004.

I. Babuska. The finite element method for elliptic equations with discontinuous
coefficients. Computing, 5:207—213, 1970.

I. Babuska and S.A. Sauter. Is the pollution effect of the fem avoidable for the
helmholtz equation considering high wave number? SIAM J. Numer. Analy.,
34:2392—2423. 1997.

F. Baetke, H. Werner, and H. Wengle. Numerical-simulation of turbulent-flow
over surface-mounted obstacles with sharp edges and corners. J. Wind Engng.
and Indust. Aerodyn, 35:129—147, 1990.

N.A. Baker. Improving implicit solvent simulations: a poisson-centric view.
Current Opin. Struct. Biol., 15:137-143, 2005.

N.A. Baker, M. Holst, and F. Wang. Adaptive multilevel finite element solution
of the poisson-boltzmann equation ii. refinement at solvent-accessible surfaces
in biomolecular systems. J. Comp. Chem, 21:1343—1352, 2000.

N.A. Baker, D. Sept, M.J. Holst, and J .A. McCammon. The adaptive multilevel
finite element solution of the poisson-boltzmann equation on massively parallel
computers. IBM J. Res. Dev., 45:427—438, 2001.

N.A. Baker, D. Sept, M.J. Holst S. Joseph, and J.A. McCammon. Electro-
statics of nanosystems: application to microtubules and the ribosome. PNAS,
98:10037—10041, 2001.

G. Bao, G.W. Wei, and S. Zhao. Numerical solution of the helmholtz equation
with high wave numbers. Int. J. Numer. Meth. Engng., 59:389—408, 2004.

P. Bates, G.W. Wei, and S. Zhao. Minimal molecular surfaces and their appli-
cations. J. Comput. Chem, 29:380—391, 2008.

RA. Berthelsen. A decomposed immersed interface method for variable co-
efficient elliptic equations with non-smooth and discontinuous solutions. J.
Comput. Phys, 197:364—386, 2004.

G. Biros, L.X. Ying, and D. Zorin. A fast solver for the stokes equations with
distributed forces in complex geometries. J. Comput. Phys, 193:317—348, 2004.

123

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24l

[25]

[26]

[27]

A. Bordner and G. Huber. Boundary element solution of the linear poisson—
boltzmann equation and a multipole method for the rapid calculation of forces
on macromolecules in solution. J. Comput. Chem, 24:353—367, 2003.

A. Boschitsch and M. Fenley. Hybrid boundary element and finite difference
method for solving the nonlinear poisson—boltzmann equation. J. Comput.
Chem, 25:935-955, 2004.

A. Boschitsch, M. Fenley, and H.-X. Zhou. Fast boundary element method
for the linear poisson—boltzmann equation. J. Phys. Chem, B 106:2741—2754,
2002.

J. Bramble and J. King. A finite element method for interface problems in
domains with smooth boundaries and interfaces. Adv. Comput. Math, 6:109—
138, 1996.

BR. Brooks, R.E. Bruccoleri, B. D. Olafson, D.J. States, S. Swaminathan, and
M. Karplus. Charmm: A program for macromolecular energy, minimization,
and dynamics calculations. J. Comput. Chem, 4:187—217, 1983.

J.F. Huang B.Z. Lu, X.L. Cheng and J.A. McCammon. Order 11 algorithm
for computation of electrostatic interactions in biomolecular systems. PNAS,
103:19314—19319, 2006.

W. Cai and S.Z. Deng. An upwinding embedded boundary method for maxwell’s
equations in media with material interfaces: 2d case. J. Comput. Phys,
190:159—183, 2003.

S. Caorsi, M. Pastorino, and M. Raffetto. Electromagnetic scattering by a con-
ducting strip with a multilayer elliptic dielectric coating. IEEE T. Electromag.
Compat., 41:335—343, 1999.

Z.J. Cendes. D.N. Shenton, and H. Shahnasser. Magnetic field computation
using delaney triangulation and complementary finite element methods. IEEE
T. Magn., 19:2251—2554, 1983.

D. L. Chapman. A contribution to the theory of electrocapollarity. Phil. Mag.,
25:475—481, 1913.

1. Chem, J .G. Liu, and WC. Wang. Electrostatics calculations: latest method-
ological advances. PMethods and Application of Analysis, 10:309—328, 2003.

ML. Connolly. Molecular surface triangulation. J. Appl. Crystallogr., 18:499—
505, 1985.

GM. Cortis and RA. Friesner. Numerical solution of the poisson-boltzmann
equation using tetrahedral finite-element meshes. J. Comput. Chem, 1811591—
1608, 1997.

ME. Davis, J .D. Madura, B.A. Luty, and J .A. McCammon. Electrostatics and
diffusion of molecules in solution: simulations with the university of houston
brownian dynamics program. Comput. Phys. Commun., 62:187— 197, 1991.

124

[28] S.Z. Deng, K. Ito, and Z.L. Li. Three-dimensional elliptic solvers for interface
problems and applications. J. Comput. Phys, 184:215—243, 2003.

[29] MA. Dumett and JP. Keener. An immersed interface method for solving
anisotropic elliptic boundary value problems in three dimensions. SIAM J. Sci.
Comput, 25:348—367, 2003.

[30] F. Eisenhaber and P. Argos. Improved strategy in analytic surface calculation
for molecular systems: Handling of singularities and computational efficiency.
J. Comput. Chem, 141272—1280, 1993.

[31] EA. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof. Combined immersed-
boundary finite-difference methods for three-dimensional complex flow simula—
tions. J. Comput. Phys, 161:30—60, 2000.

[32] RP. Fedkiw, T. Aslam, B. Merriman, and S. Osher. A non-oscillatory eule-
rian approach to interfaces in multimaterial flows (the ghost fluid method). J.
Comput. Phys, 152:457—492, 1999.

[33] M. Feig and C. L. Brooks III. Recent advances in the development and ap-
plication of implicit solvent models in biomolecule simulations. Current Opin.
Struct. Biol, 14:217—224, 2004.

[34] M. Feig, A. Onufriev, M.S. Lee, D.A. Case W. Im, and CL. Brooks III. Perfor-
mance comparison of generalized born and poisson methods in the calculation of
electrostatic solvation energies for protein structures. J. Comp. Chem, 25:265—
284, 2004.

[35] AL. Fogelson and J .P. Keener. Immersed interface methods for neumann and
related problems in two and three dimensions. SIAM J. Sci. Comput, 22:1630—
1654, 2000.

[36] F. Fogolari, A. Brigo, and H. Molinari. The poisson-boltzmann equation for
biomolecular electrostatics: a tool for structural biology. J. Mol. Recognit,
15:377—392, 2002.

[37] F. Fogolari, A. Brigo, and H. Molinari. Protocol for mm/pbsa molecular dy-
namics simulations of proteins. Biophys. J., 85:159—166, 2003.

[38] KB]. Forsten, R.E. Kozack, D.A. Lauffenburger, and S. Subramaniam. Nu-
merical solution of the nonlinear poisson-boltzmann equation for a membrane-
electrolyte system. J. Phys. Chem, 98:5580~5586, 1994.

[39] M. Francois and W. Shyy. Computations of drop dynamics with the immersed
boundary method, part 2: Drop impact and heat transfer. Numer. Heat Trans.
Part B-F’und., 44:119—143, 2003.

[40] M. Tong G. Pan and B. Gilbert. Multiwavelet based moment method under
discrete sobolev-type norm. Microwave Opt Techn. Lett., 40:47—50, 2004.

[41] W.H. Geng and G.W. Wei. Interface technique based molecular dynamics
method via the poisson-boltzmann euqtaion. J. Compat. Phys, submitted,
2008.

125

 

 

 

[42]

[431

[44]

[45]

[46]

[47]

[48]

[49]

[51]

[52]

[53]

[54]

[55]

[56]

W.H. Geng, S.N. Yu, and G.W. Wei. Treatment of charge singularities in
implicit solvent models. J. Chem. Phys, 1272114106, 2007.

RE. Georgescu, E.G. Alexov, and MR. Gunner. Combining conformational
flexibility and continuum electrostatics for calculating pkas in proteins. Biophys.
J., 83:1731—1748, 2002.

F. Gibou and R. P. Fedkiw. A fourth order accurate discretization for the
laplace and heat equations on arbitrary domains, with applications to the stefan
problem. J. Comput. Phys, 202:577—601, 2005.

M.k. Gilson, ME. Davis, B.A. Luty, and J.A. McCammon. Computation of
electrostatic forces on solvated molecules using the poisson-boltzmann equation.
J. Phys. Chem, 97:3591-—3600, 1993.

V. Gogonea and E. Osawa. Implementation of solvent effect in molecular
mechanics. 1. model development and analytical algorithm for the solvent -
accessible surface area. Supramol. Chem, 3:303—317, 1994.

G. Gouy. Constitution of the electric charge at the surface of an electrolyte. J.
Phys, 9:457—468, 1910.

J .A. Grant, B.T. Pickup, and A. Nicholls. A smooth permittivity function for
poisson-boltzmann solvation methods. Journal of Computational Chemistry,
22:608—40, 2001.

G. Guyomarch and C.-O. Lee. A discontinuous galerkin method for elliptic
interface problems with application to electroporation. AIST DAM Research
Report, page 39186, 2004.

CR. Hadley. High-accuracy finite-difference equations for dielectric waveguide
analysis i: uniform regions and dielectric interfaces. Journal of Lightwave T ech-
nology, 20:1210—1218, 2002.

J.S. Hesthaven. High-order accurate methods in time-domain computational
electromagnetics. a review. Advances in Imaging and Electron Physics, 127:59—
123, 2003.

M. Holst, N .A. Baker, and F. Wang. Adaptive multilevel finite element solution
of the poisson-boltzmann equation i. algorithms and examples. J. Comput.
Chem., 21:1319—1342, 2000.

M. Holst and F. Saied. Multigrid solution of the poisson-boltzmann equation.
J. Comput. Chem., 14:105—113, 1993.

M.J. Holst. The Poisson-Boltzmann Equation. PhD thesis, Numerical Com-
puting Group, University of Illinois, Urbana-Champaign, 1993.

B. Honig and A. Nicholls. Classical electrostatics in biology and chemstry.
Science, 268:1144-1149, 1995.

S. Hou and X.-D. Liu. A numerical method for solving variable coefficient
elliptic equation with interfaces. J. Comput. Phys, 202:411-—445, 2005.

126

 

 

[57] T.Y. Hou, Z.L. Li, S. Osher. and H. Zhao. A hybrid method for moving interface
problems with application to the hele—shaw flow. J. Comput. Phys, 134:236--
252, 1997.

[58] H. Huang and Z.L. Li. Convergence analysis of the immersed interface method.
IMA J. Numer. Anal, 19:583—608, 1999.

[59] J .K. Hunter, Z.L. Li, and H. Zhao. Reactive autophobic spreading of drops. J.
Comput. Phys, 183:335—366, 2002.

[60] G. Iaccarino and R. Verzicco. Immersed boundary technique for turbulent flow
simulations. Appl. Mech. Rea, 56:331—347, 2003.

[61] W. Im, D. Beglov, and B. Roux. Continuum solvation model: Computation of
electrostatic forces from numerical solutions to the poisson-boltzmann equation.
Comput. Phys. Commun., 111:59—75, 1998.

[62] S. Jin and X.L. Wang. Robust numerical simulation of porosity evolution in

chemical vapor infiltration ii. two—dimensional anisotropic fronts. J. Comput.
Phys, 179:557—577, 2002.

[63] H. Johansen and P. Colella. A cartesian grid embedding boundary method for
poisson’s equation on irregular domains. J. Comput. Phys, 147:60—85, 1998.

[64] A. D. MacKerell jr., D. Bashford, M. Bellott, J. D. Dunbrack, M. J. Evanseck,
M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir,
K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen,
B. Prodhom, and W. E. Reiher. All-atom empirical potential for molecular
modeling and dynamics studies of proteins. J. Phys. Chem., 102:3586—3616,
1998.

[65] A. Juffer, E. Botta, B. van Keulen, A. van der Ploeg, , and H. Berendsen. The
electric potential of a macromolecule in a solvent: a fundamental approach. J.
Comput. Phys, 97:144—171, 1991.

[66] J .D. Kandilarov. Immersed interface method for a reaction-diffusion equation
with a moving own concentrated source. Lecture Notes Comput. Sci, 2542:506—
513, 2003.

[67] J .G. Kirkwood. Theory of solution of molecules containing widely separated
charges with special application to zwitterions. J. Chem. Phys, 7:351—361,
1934.

[68] I. Klapper, R. Hagstrom, R. Fine, K. Sharp, and B. Honig. Focusing of electric
fields in the active site of cu-zn superoxide dismutase: Effects of ionic strength
and amino-acid modification. Proteins, 1:47—59, 1986.

[69] P. Koehl. Electrostatics calculations: latest methodological advances. Current
Opin. Struct. Biol, 16:142—151, 2006.

[70] PA. Kollman, I. Massova, C. Reyes, B. Kuhn, S. Huo, L. Chong, M. Lee,
T. Lee, Y. Duan, and W. Wang et al. Calculating structuras and free energies
of complex molecules: combining molecular mechanics and continuum models.
Acc. Chem. Res, 33:889—897, 2000 .

127

 

 

[71] B. Lee and F .M. Richards. Interpretation of protein structures - estimation of
static accessibility. J. Mol Biol, 55:379, 1971.

[72] L. Lee and R.J. LeVeque. An immersed interface method for incompressible
navier-stokes equations. SIAM J. Sci. Comput., 25:832—856, 2003.

[73] R.J. LeVeque and Z.L. Li. The immersed interface method for elliptic equations
with discontinuous coefficients and singular sources. SIAM J. Numer. Anal,
31:1019—1044, 1994.

[74] Z.L. Li. A fast iterative algorithm for elliptic interface problems. SIAM J.
Numer. Anal, 35:230—254, 1998.

[75] Z.L. Li and K. Ito. Maximum principle preserving schemes for interface prob-
lems with discontinuous coefficients. SIAM J. Sci. Comput., 23:339—361, 2001.

[76] Z.L. Li, T. Lin, and X.H. Wu. New cartesian grid methods for interface problems
using the finite element formulation. Numerische Mathematik, 96:61—98, 2003.

[77] Z.L. Li and SR. Lubkin. Numerical analysis of interfacial two—dimensional
stokes flow with discontinuous viscosity and variable surface tension. Int. J.
Numer. Meth. F laid, 37:525—540, 2001.

[78] Z.L. Li, W.-C. Wang, I.-L. Chem, and M.-C. Lai. New formulations for interface
problems in polar coordinates. SIAM J. Sci. Comput., 25:224—245, 2003.

[79] J. Liang and S. Subranmaniam. Computation of molecular electrostatics with
boundary element methods. Biophst., 73:1830—1841, 1997.

[80] M. N. Linnick and H. F. Fasel. A high-order immersed interface method for
simulating unsteady incompressible flows on irregular domains. J. Comput.
Phys, 204:157—192, 2005.

[81] W.K. Liu, D. Farrell Y. Liu, X. Wang L. Zhang, N. Patankar Y. Fukui,
Y. Zhang, C. Bajaj, X. Chen, and H. Hsu. Immersed finite element method and
its applications to biological systems. Computer Methods in Applied Mechanics
and Engineering, 195:1722—1749, 2006.

[82] X.D. Liu, R.P. Fedkiw, and M. Kang. A boundary condition capturing method
for Poisson’s equation on irregular domains, page J. Comput. Phys, 160.

[83] B. Lombard and J. Piraux. How to incorporate the spring-mass conditions in
finite—difference schemes. SIAM J. Sci. Comput., 24:1379—1407, 2003.

[84] B.Z. Lu, W.Z. Chen, C.X. Wang, and X.J. Xu. Protein molecular dynam-
ics with electrostatic force entirely determined by a single poisson-boltzmann
calculation. Proteins, 48:497—504, 2002.

[85] Q. Lu and R. Luo. A poisson-boltzmann dynamics method with nonperiodic
boundary condition. J. Chem. Phys, 119:11035—11047, 2003.

[86] R. Luo, L. David, and MK. Gilson. Accelerated poisson-boltzmann calculations
for static and dynamic systems. J. Comput. Chem., 23:1244—1253, 2002.

128

 

 

[87]

[88]

[89]

[90]

[91]

[92]

[93]

[94]

[95]

[96]

[97]

[98]

[99]

[100]

[101]

BB. Macak, W.D. Munz, and J .M. Rodenburg. Plasma-surface interaction at
sharp edges and corners during ion-assisted physical vapor deposition. part i:

Edge-related effects and their influence on coating morphology and composition.
J. A ppl. Phys, 94:2829-2836, 2003.

A. Mayo. The fast solution of poisson’s and the biharmonic equations on irreg-
ular regions. SIAM J. Numer. Anal, 21:285—299, 1984.

A. Mckenney, L. Greengard, and A. Mayo. A fast poisson solver for complex
geometries. J. Comput. Phys, 118:348—355, 1995.

R. Miniowitz and J .P. Webb. Covariant-projection quadrilateral elements for
the analysis of wave-guides with sharp edges. IEEE T. Microwave Theory and
Techn., 39:501—505, 1991.

R. Mittal and G. Iaccarino. Immersed boundary methods. Annu. Rev. Fluid
Mech., 37:236—261, 2005.

J.E. Nielsen and J.A. McCammon. On the evaluation and optimization of
protein x-ray structures for pka calculations. Protein Sci, 12:313—326, 2003.

Z. Pantic-Tanner, J.Z. Savage, D.R. Tanner, and AF. Peterson. Two-
dimensional singular vector elements for finite-element analysis. IEEE T. Mi-
crowave Theory Techn., 46:178—184, 1998.

CS. Peskin. Numerical analysis of blood flow in heart. J. Comput. Phys,
25:220—252, 1977.

CS. Peskin. Lectures on mathematical aspects of physiology. Lectures in Appl.
Math, 19269—107, 1981.

CS. Peskin and D.M. Mcqueen. A 3-dimensional computational method for
blood-flow in the heart. 1.immersed elastic fibers in a viscous incompressible
fluid. J. Comput. Phys, 81:372—405, 1989.

N .V. Prabhu, P. Zhu, and K.A. Sharp. Implementation and testing of stable,
fast implicit solvation in molecular dynamics using the smooth-permittivity
finite difference poisson-boltzmann method. J. Comput. Chem., 25:2049—2064,
2004.

Z. H. Qiao, Z.L. Li, and T. Tang. A finite difference scheme for solving the
nonlinear poisson-boltzmann equation modeling charged spheres. J. Comput.
Chem., 24:252—264, 2006.

F. M. Richards. Areas, volumes, packing and protein structure. Annu. Rev.
Biophys. Bioeng, 6:151-176, 1977.

W. Rocchia, E. Alexov, and B. Honig. Extending the applicability of the nonlin-
ear poisson-boltzmann equation: multiple dielectric constants and multivalent
ions. J.Phys. Chem. B, 105:6507—6514, 2001.

B. Roux and T. Simonson. Implicit solvent models. Biophys. Chem., 78:39102,
1999.

129

 

 

[102] M.F. Sanner, A.J. Olson, and J .C. Spehner. Reduced surface: An efficient way
to compute molecular surfaces. Biopolymers, 38:305—320, 1996.

[103] M. Schulz and G. Steinebach. Two-dimensional modelling of the river rhine. J.
Comput. Appl. Math, 145:39406, 2002.

[104] J .A. Sethian and A. Wiegmann. Structural boundary design via level set and
immersed interface methods. J. Comput. Phys, 163:489—528, 2000.

[105] K.A. Sharp and B. Honig. Electrostatic interactions in macromolecules: thoery
and applications. Annu. Rev. Biophys. Biophys. Chem., 19:301—332, 1990.

[106] P. Smereka. The numerical approximation of a delta function with application
to level set methods. J. Comput. Phys, 211:77—90, 2006.

[107] J .M.J . Swanson, R.H. Henchman, and J .A. McCammon. Revisiting free energy
calculations: a theoretical connection to mm/pbsa and direct calculation of the
association free energy. Biophys. J., 86:67-74, 2004.

[108] J.M.J. Swanson, J. Mongan, and J.A. McCammon. Limitations of atom-
centered dielectric functions in implicit solvent models. J. Phys. Chem.,
109:14769—14772, 2005.

[109] J.M.J. Swanson, J.A. Wagoner, N.A. Baker, and J.A. McCammon. Optimiz-
ing the poisson dielectric boundary with explicit solvent forces and energies:
lessons learned with atom-centered dielectric functions. J Chem Theory Com-
put., 3:170—183, 2007.

[110] A.K. Tornberg and B. Engquist. Numerical approximations of singular source
terms in differential equations. J. Comput. Phys, 200:462-488, 2004.

[111] B.J.E. van Rens, W.A.M. Brekelmans, and F .P.T. Baaijens. Modelling friction
near sharp edges using a eulerian reference frame: application to aluminum
extrusion. Int. J. Numer. Methods Engng., 54:453—471, 2002.

[112] CL. Vizcarra and S.L. Mayo. Electrostatics in computational protein design.
Current Opin. Struct. Biol, 9:622—626, 2005.

[113] J .V. Voorde, J. Vierendeels, and E. Dick. Flow simulations in rotary volumetric
pumps and compressors with the fictitious domain method. J. Comput. Appl.
Math, 168:491—499, 2004.

[114] Y.N. Vorobjev and H .A. Scheraga. A fast adaptive multigrid boundary element
method for macromolecular electrostatic computations in a solvent. J. Comput.
Chem., 18:569—583, 1997.

[115] J. Wagoner and N .A. Baker. Solvation forces on biomolecular structures: A
comparison of explicit solvent and poisson-boltzmann models. J. Comput.
Chem., 25:1623—1629, 2004.

[116] J .H. Walther and G. Morgenthal. An immersed interface method for the vortex-
in—cell algorithm. J. Turbulence, 3:Art. No. 039, 2002.

130

 

 

[117] J. Warwicker. Improved pka calculations through flexibility based sampling of
a water-dominated interaction scheme. Protein Sci., 13:2793—2805, 2004.

[118] J. Warwicker and HG. Watson. Calculation of the electric-potential in the
active-site cleft due to alpha-helix dipoles. J. Mol Biol, 154:671—679, 1982.

[119] A. Wiegmann and KP. Bube. The explicit-jump immersed interface method:
finite difference methods for pdes with piecewise smooth solutions. SIAM J.
Numer. Anal, 37:827-862, 2000.

[120] SN. Yu. Matched Interface and Boundary (MIB) Method for Geometric Singu-
larities and Its Application to Molecular Biology And Structural Analysis. PhD
thesis, Department of Mathematics, Michigan State University, 2007.

[121] SN Yu, W.H. Geng, and G. W. Wei. Treatment of geometric singularities in
the implicit solvent models. J. Chem. Phys, 126:244108, 2007.

[122] SN. Yu and G.W. Wei. Three-dimansional matched interface and boundary
(MIB) method for treating geometric singularities. J. Comput. Phys, submit-
ted, 2007.

[123] SN. Yu, Y.C. Zhou, and G.W. Wei. Matched interface and boundary (MIB)
method for elliptic problems with sharp-edged interfaces. J. Comput. Phys,
224:729—756, 2007.

[124] R.J. Zauhar and RS. Morgan. A new method for computing the macromolec—
ular electric-potential. J. Mol. Biol, 186:815—820, 1985.

[125] Y. Zhang, G. Xu, and C. Bajaj, Quality meshing of implicit solvation models of
biomolecular structures. Computer Aided Geometric Design, 23:510—530, 2006.

[126] S. Zhao and G.W. Wei. High order F DTD methods via derivative matching for
maxwell’s equations with material interfaces. J. Comput. Phys, 200:60—103,
2004.

[127] Y.C. Zhou. Matched Interface and Boundary(MIB) Method and its Applications
to the implicit solvent model. PhD thesis, Department of Mathematics, Michigan
State University, 2006.

[128] Y.C. Zhou, M. Feig, and G.W. Wei. Highly accurate biomolecular electrostatics
in continuum dielectric environments. J. Comput. Chem., 29:87—97, 2008.

[129] Y.C. Zhou and G.W. Wei. On the fictitious-domain and interpolation formula-
tions of the matched interface and boundary (MIB) method. J. Comput. Phys,
219:228—246, 2006.

[130] Y.C. Zhou, S. Zhao, M. Feig, and G. W. Wei. High order matched interface and
boundary (MIB) schemes for elliptic equations with discontinuous coefficients
and singular sources. J. Comput. Phys, 213:39112, 2006.

[131] Z. Zhou, P. Payne, M. Vasquez, N. Kuhn, and M. Levitt. F mite-difference so-
lution of the poisson-boltzmann equation: Complete elimination of self—energy.
J. Comput. Chem., 17:1344-1351, 1996.

131

 

[ll]

58

       

ur]l]l]]]]]]]]i[i]l]i]