Corrected Trapezoidal Rule-IBIM for linearized Poisson-Boltzmann equation

2025-05-06 0 0 1.97MB 22 页 10玖币

侵权投诉

Corrected Trapezoidal Rule-IBIM for linearized

Poisson-Boltzmann equation

Federico Izzo∗,†Yimin Zhong‡Olof Runborg§Richard Tsai¶

Abstract

In this paper, we solve the linearized Poisson-Boltzmann equation, used to model the

electric potential of macromolecules in a solvent. We derive a corrected trapezoidal rule with

improved accuracy for a boundary integral formulation of the linearized Poisson-Boltzmann

equation. More speciﬁcally, in contrast to the typical boundary integral formulations, the

corrected trapezoidal rule is applied to integrate a system of compacted supported singular

integrals using uniform Cartesian grids in R3, without explicit surface parameterization. A

Krylov method, accelerated by a fast multipole method, is used to invert the resulting linear

system. We study the eﬃcacy of the proposed method, and compare it to an existing, lower

order method. We then apply the method to the computation of electrostatic potential of

macromolecules immersed in solvent. The solvent excluded surfaces, deﬁned by a common

approach, are merely piecewise smooth, and we study the eﬀectiveness of the method for such

surfaces.

Key words: Poisson-Boltzmann equation; implicit boundary integral method; implicit

solvent model; singular integrals; trapezoidal rules.

AMS subject classiﬁcations 2020: 45A05, 65R20, 6N5D30, 65N80, 78M16, 92E10

1 Introduction

We propose to solve the linearized Poisson-Boltzmann equation using an Implicit Boundary Inte-

gral Method (IBIM), discretized by a corrected trapezoidal rule (CTR). The linearized Poisson-

Boltzmann equation can be used to model the electrostatics of charged macromolecule-solvent

systems. Numerical modeling and simulation of such systems is an active and relevant research

topic. Electrostatic interactions with the solvent signiﬁcantly aﬀect the overall macromolecule be-

havior. They are relevant for example in electrochemistry, as an aqueous solvent plays a signiﬁcant

role in the dynamical processes of biological molecules (see [1,4] for comprehensive introductions

to the problem).

When describing the electrostatic interactions, implicit solvent methods approximate the prob-

lem by simplifying the description of the solvent environment and consequently greatly reducing

the degrees of freedom. The frameworks for approximating implicitly the electrostatic interactions

are for example the Coulomb-ﬁeld approximation, the generalized Born models, and the Poisson-

Boltzmann theory. The latter two are more popular, and among them the Poisson-Boltzmann

theory oﬀers a more detailed representation at the cost of more expensive computations (see also

[3]).

The Poisson-Boltzmann equation, which describes the electrostatic ﬁeld determined by the

interactions in the Poisson-Boltzmann theory, is nonlinear. The solution to its linearized form is

however widely used, both because it represents a valid approximation to the nonlinear solution

in certain settings, and because it can be used iteratively to ﬁnd the original nonlinear solution.

∗Corresponding author

†Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden (izzo@kth.se)

‡Department of Mathematics and Statistics, Auburn University, Auburn AL, USA (yimin.zhong@auburn.edu)

§Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden (olofr@kth.se)

¶Department of Mathematics and Oden Institute for Computational Engineering and Sciences, The University

of Texas at Austin, Austin TX, USA (ytsai@math.utexas.edu)

arXiv:2210.03699v1 [math.NA] 7 Oct 2022

Hence, numerically solving the linearized Poisson-Boltzmann equation is still an active and relevant

problem (see e.g. [11,7,18,23,22]).

In this work, we propose a fast numerical method for solving the linearized Poisson-Boltzmann

equation on surfaces:

−∇ · (εI∇ψ(x)) = PNc

k=1 qkδ(x−zk),in Ω,

−∇ · (εE∇ψ(x)) = −¯κ2

Tψ(x) in ¯

Ωc,

ψ(x)|Γ−=ψ(x)|Γ+on Γ,

εI

∂ψ

∂n(x)Γ−

=εE

∂ψ

∂n(x)Γ+

on Γ,

|x|ψ(x)→0,|x|2|∇ψ(x)| → 0,as |x|→∞.

(1.1)

In (1.1), ψrepresents the electrostatic potential, Ncis the number of atoms composing the macro-

molecules, which have centers {zk}Nc

k=1, radii {rk}Nc

k=1, and charge numbers {qk}Nc

k=1 respectively.

Γ represents the closed surface which separates the region occupied by the macromolecule (here

denoted by Ω) and the rest of the space ¯

Ωc:= R3\¯

Ω = R3\(Ω∪Γ). We will describe in Section 2.6

how we test the method on simple smooth surfaces Γ, and then in Section 3how we apply it to

actual solvent-molecule interfaces. The operator ∂ψ

∂n(x) represents the normal derivative of ψin

x∈Γ with normal pointing outward from Γ. The parameters εIand εEare the dielectric constants

inside and outside Ω respectively, and ¯κTis a screening parameter. The radiation conditions on

the last row are needed to ensure uniqueness of the solutions.

We use the following boundary integral formulation [11] to ﬁnd the solution to the equations

(1.1):

21 + εE

εIψ(x) + ZΓ

K11(x,y)ψ(y)dy−ZΓ

K12(x,y)∂ψ

∂n(y)dy=

k=1

εI

G0(x,zk),(1.2)

21 + εI

εE∂ψ

∂n(x) + ZΓ

K21(x,y)ψ(y)dy−ZΓ

K22(x,y)∂ψ

∂n(y)dy=

k=1

εI

∂G0

∂nx

(x,zk).(1.3)

where

K11(x,y) := ∂G0

∂ny

(x,y)−εE

εI

∂Gκ

∂ny

(x,y), K12(x,y) := G0(x,y)−Gκ(x,y),

K21(x,y) := ∂2G0

∂nx∂ny

(x,y)−∂2Gκ

∂nx∂ny

(x,y), K22(x,y) := ∂G0

∂nx

(x,y)−εI

εE

∂Gκ

∂nx

(x,y).

(1.4)

The fundamental solutions which deﬁne the kernels in (1.4) are

G0(x,y) := 1

4π|x−y|, Gκ(x,y) := e−κ|x−y|

4π|x−y|.(1.5)

The integrals in (1.2)-(1.3) are well deﬁned because of the integrability of the kernels for x=y:

K11(x,y)∼ O(|x−y|−1), K21(x,y)∼ O(|x−y|−1), K22(x,y)∼ O(|x−y|−1),(1.6)

K12(x,y)∼ O(1).

As in [22], we write these integrals in the Implicit Boundary Integral Method (IBIM) framework

(see [13,12]).

We point out that the kernel K21 is integrable even though it is deﬁned as the diﬀerence between

two hypersingular kernels which are not Cauchy integrable. This is the power of the boundary

integral formulation derived by Juﬀer et al. [11]. The kernel K21 proves to be the bottleneck of

the accuracy order attainable given a second order approximation of the surface around the target

point x. If the approximation is improved to third order, it is possible to improve the order of

accuracy for the kernels K11,K22, and K12 but not for K21. Only a fourth order approximation

allows us to increase the order of accuracy for K21. This will be explained in more detail in

Section 2.4.1.

2 Numerical solution of the volumetric boundary integral

formulation

In this section we will provide the necessary background information for the numerical solution of

(1.2)-(1.3).

In Section 2.1 we present the volumetric formulation of the surface integrals via non-parametrized

surface representation. In Section 2.2 we present the discretization (2.14)-(2.15) of the original BIEs

(1.2)-(1.3) using the framework shown in the previous subsection. In Section 2.3 we brieﬂy present

the kernel regularization (K-reg) developed in [22]. In Section 2.4 we present the formally second

order accurate quadrature rule (CTR2) based on the trapezoidal rule which we use to approximate

the singular volume integrals. In Section 2.5 we present an overview of the algorithms used to

apply the quadrature rule given the surface Γ. In Section 2.6 we present numerical tests on smooth

surfaces to show the order of accuracy of CTR2 and compare it to K-reg.

2.1 Extensions of the singular boundary integral operators

Let Ω ⊂R3be a bounded open set with C2boundaries, and ∂Ω =: Γ. We shall refer to Γ as the

surface. Let fbe a function deﬁned on Γ. In this section we present an approach for extending a

boundary integral ZΓ

f(y)dσy,(2.1)

to a volumetric integral around the surface. Instead of parameterizations, this approach relies on

the Euclidean distance to the surface, and its derivatives. More precisely, we deﬁne the signed

distance function

dΓ(x) := (miny∈Γkx−yk,if x∈Ω

−miny∈Γkx−yk,if x∈Ωc(2.2)

and the closest point projection

PΓ(x) := argminy∈Γkx−yk.(2.3)

If there is more than one global minimum, we pick one randomly from the set. Let CΓdenote

the set of points in R3which are equidistant to at least two distinct points on Γ. The reach τis

deﬁned as infx∈Γ,y∈CΓkx−yk. Clearly, τis restricted by the local geometry (the curvatures) and

the global structure of Γ (the Euclidean and geodesic distances between any two points on Γ).

In this paper, we assume that Γ is C2and thus has a non-zero reach. Let Tεdenote the set of

points of distance at most εfrom Γ:

Tε={x∈R3:|dΓ(x)| ≤ ε}.(2.4)

Then PΓis a diﬀeomorphism between Γ and the level sets of dΓin T. Furthermore,

PΓ(x) = x−dΓ(x)∇dΓ(x),x∈Tε.

We deﬁne the extension of the integrand fby

f(x) := f(PΓx),x∈R3.(2.5)

As in [12,13], we can then rewrite the surface integral (2.1):

ZΓ

f(y)dσy=ZR3

f(x)δΓ,ε(y)dy=ZTε

f(x)δΓ,ε(y)dy,(2.6)

where

δΓ,ε(y) := JΓ(y)δε(dΓ(y)),y∈R3,

δε(η) := 1

εφη

ε, φ ∈C∞(R) supported in [−1,1],and ZR

φ(x)dx= 1,(2.7)

and JΓ(y0) is the Jacobian of the transformation from Γ to the level set Γη:= {y∈Rn:dΓ(y) =

η}.In R3,JΓ(y0) is a quadratic polynomial in dΓ(y0):

JΓ(y0) := 1 + 2dΓ(y0)H(y0) + dΓ(y0)2G(y0).(2.8)

where H(y0) and G(y0) are respectively the mean and Gaussian curvatures of Γηat y0. The

curvatures as well as the corresponding principal directions can be extracted from PΓ; see [13] for

more detail.

Our primary focus is when f(y) is replaced by a function of the kind K(x,y)ζ(y), corresponding

to the kernels and potentials found in (1.2)-(1.3), (1.4)), and the integral operators of form

Ji[ζ1, ζ2](x) := ZΓ

Ki1(x,y)ζ1(y)dy−ZΓ

Ki2(x,y)ζ2(y)dy,x∈R3, i = 1,2.(2.9)

We can now write the operators (2.9) in volumetric form:

Ji[ρ1, ρ2](x) := ZTε

Ki1(x,y)ρ1(y)δΓ,ε(y)dy−ZTε

Ki2(x,y)ρ2(y)δΓ,ε(y)dy,x∈T, i = 1,2,

(2.10)

where

K(x,y) := K(x, PΓy),x,y∈Rn.(2.11)

It is important to notice that if K(x,y) is singular for x=y, then K(PΓx,y) is singular on the

set

{(x,y)∈Rn×Rn:PΓx=PΓy},

i.e. for a ﬁxed x∗∈Γ, K(x∗,y) is singular along the normal line passing through x∗, while

K(x∗,y) is singular in a point.

Finally, with

λ1:= 1

21 + εE

εI, λ2:= 1

21 + εI

εE,

g1(x) := PNc

k=1

εI

G0(x,zk), g2(x) := PNc

k=1

εI

∂G0

∂nx

(x,zk),

(2.12)

the volumetric forms of equations (1.2)-(1.3) are derived:

λiρi(x) + Ji[ρ1, ρ2](PΓx) = gi(PΓx),x∈Tε, i = 1,2.(2.13)

The solutions ρ1,ρ2will coincide with the constant extensions along the normals of ψand ψn

respectively: ρ1(x)≡ψ(PΓ(x)), ρ2(x)≡ψn(PΓ(x)); see, for example, the arguments given in [9].

2.2 Quadrature rules on uniform Cartesian grids

In this paper, we derive numerical quadratures for the singular integral operators Ji[ρ1, ρ2](x)

for x∈Γ (equivalently, Ji[ρ1, ρ2](PΓx) for x∈Tε) for i= 1,2. The quadrature rules will be

constructed based on the trapezoidal rule for the grid nodes Th

ε:= Tε∩hZ3, which corresponds to

the portion of the uniform Cartesian grid hZ3within Tε. Since the integrand in (2.10) is singular

for x∈Γ, the trapezoidal rule should be corrected near xfor faster convergence. Correction will

be deﬁned by summing the judiciously derived weights over a set of grid nodes denoted by Nh(x).

The sum will be denoted by Rh(x). Ultimately, the quadrature for Ji[ρ1, ρ2](PΓx) will involve

the regular Riemann sum of the integrand in Th

ε\Nh(x), and the correction Rh(x) in Nh(x):

J1[ρ1, ρ2](x)≈h3X

ym∈Th

ε\Nh(x)

K11(x,ym)ρ1(ym)δΓ,ε(ym) + h2X

ym∈Nh(x)

ω(11)

mρ1(ym)δΓ,ε(ym)

−h3X

ym∈Th

ε\Nh(x)

K12(x,ym)ρ2(ym)δΓ,ε(ym)−h3X

ym∈Nh(x)

ω(12)

mρ2(ym)δΓ,ε(ym),

and

J2[ρ1, ρ2](x)≈h3X

ym∈Th

ε\Nh(x)

K21(x,ym)ρ1(ym)δΓ,ε(ym) + h2X

ym∈Nh(x)

ω(21)

mρ1(ym)δΓ,ε(ym)

−h3X

ym∈Th

ε\Nh(x)

K22(x,ym)ρ2(ym)δΓ,ε(ym)−h2X

ym∈Nh(x)

ω(22)

mρ2(ym)δΓ,ε(ym),

where ω(ij)

mis a weight which depends on the kernel Kij , on the grid node ym, and on the principal

curvatures and directions of the surface in x∈Γ.

We now have all the ingredients to write out the linear system we need to solve: ∀yk∈Th

ε,

λ1ρ1(yk)+h3X

ym∈Th

ε\Nh(yk)K11(PΓyk,ym)ρ1(ym)−K12(PΓyk,ym)ρ2(ym)δΓ,ε(ym)

+h2X

ym∈Nh(yk)ω(11)

mk ρ1(ym)−h ω(12)

mk ρ2(ym)δΓ,ε(ym) = g1(PΓyk),(2.14)

λ2ρ2(yk)+h3X

ym∈Th

ε\Nh(yk)K21(PΓyk,ym)ρ1(ym)−K22(PΓyk,ym)ρ2(ym)δΓ,ε(ym)

+h2X

ym∈Nh(yk)ω(21)

mk ρ1(ym)−ω(22)

mk ρ2(ym)δΓ,ε(ym) = g2(PΓyk),(2.15)

where ρ1(yk)≈¯

ψ(yk) and ρ2(yk)≈¯

ψn(yk). Corresponding to the target node ykthe set Nh(yk)

contains all the correction nodes, and ω(ij)

mk =ω[s(ij);αmk, βmk], i, j = 1,2, are the correction

weights, dependent on the kernel they correct1via the function s(ij), and on (αmk, βmk), parameters

dependent on h,ym, and yk. The details of how the corrections nodes are chosen are presented in

Section 2.4, while the deﬁnition and computation of the weights is presented in Section 2.4.2.

The contribution of this paper are two quadrature rules. One is a high order, trapezoidal

rule-based, quadrature rule for Jvia J, eﬀective for smooth surfaces which are globally C2. The

second is a hybrid rule which combines the previous trapezoidal rule-based with the constant

regularization from [22].

2.3 Kernels regularization in IBIM - K-reg

The principle of the regularization proposed in [22] is to substitute the kernel Kwith a regularized

version ˜

Karound the singularity point, so that the integral over the domain would be as close as

possible to the original one. Thus the technique is independent of the quadrature rule. Given a

parameter τ > 0 and x,y∈Tε, the regularization of the kernel is

Kτ(x,y) := (K(x,y),if |x−y|P≥τ,

CΓ,τ ,if |x−y|P< τ, (2.16)

where |x−y|Pis the distance between the projections of xand yon the tangent plane at PΓx,

and CΓ,τ is a constant dependent on the surface Γ and on the parameter τ. The constant CΓ,τ

is constructed so that it behaves in the same way as Kin a neighborhood of xdependent on τ.

Speciﬁcally, given V(PΓx;τ) disc of radius τon the tangent plane of Γ at PΓx, it is deﬁned as

CΓ,τ =1

V(PΓx;τ)ZV(PΓx;τ)

K(x, PΓz)dσz.

Then, for the kernels (1.4), the constants are

C(11)

Γ,τ =C(22)

Γ,τ =C(21)

Γ,τ = 0, C(12)

Γ,τ =exp(−κτ)−1 + κτ

2πκτ2,

for K11,K22,K21, and K12 respectively.

1Please note that the weight for the kernel K12 has a diﬀerent scaling compared to the other three. This is

because in order to reach order of accuracy two for such kernel (see (1.6)) the correction weight is zero, so we apply

the correction necessary to reach order three.

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

10 玖币 0人已下载

立即下载

摘要：

CorrectedTrapezoidalRule-IBIMforlinearizedPoisson-BoltzmannequationFedericoIzzo*;YiminZhongOlofRunborg§RichardTsai¶AbstractInthispaper,wesolvethelinearizedPoisson-Boltzmannequation,usedtomodeltheelectricpotentialofmacromoleculesinasolvent.Wederiveacorrectedtrapezoidalrulewithimprovedaccuracyforabo...

展开>> 收起<<

Corrected Trapezoidal Rule-IBIM for linearized Poisson-Boltzmann equation.pdf

共22页,预览5页

还剩页未读，继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

Corrected Trapezoidal Rule-IBIM for linearized Poisson-Boltzmann equation

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: