A Nonlinear Elasticity Model of Macromolecular Conformational Change Induced by Electrostatic Forces

0.55 MB·English

by  Yongcheng Zhou

A NONLINEAR ELASTICITY MODEL OF MACROMOLECULAR
CONFORMATIONAL CHANGE INDUCED BY ELECTROSTATIC FORCES
YONGCHENGZHOU,MICHAELHOLST,ANDJAMESANDREWMCCAMMON
ABSTRACT. In this paper we propose a nonlinear elasticity model of macromolecular
conformational change (deformation) induced by electrostatic forces generated by an
implicit solvation model. The Poisson-Boltzmann equation for the electrostatic poten-
tialisanalyzedinadomainvaryingwiththeelasticdeformationofmolecules,andanew
0 continuousmodeloftheelectrostaticforcesisdevelopedtoensuresolvabilityofthenon-
1 linearelasticityequations. Wederivetheestimatesofelectrostaticforcescorresponding
0 to four types of perturbations to an electrostatic potential field, and establish the exis-
2
tanceofanequilibriumconfigurationusingafixed-pointargument,undertheassumption

n thatthechangeintheionicstrengthandchargesduetotheadditionalmoleculescausing
a
thedeformationaresufficientlysmall. Theresultsarevalidforelasticmodelswitharbi-
J
trarilycomplexdielectricinterfacesandcavities,andcanbegeneralizedtolargeelastic
8
deformationcausedbyhighionicstrength, largecharges, andstrongexternalfieldsby

] usingcontinuationmethods.
P
A
.
h
t
a CONTENTS
m
1. AnElectro-ElasticModelofConformationalChange 2
[
2. NotationandSomeBasicEstimates 5

1 3. NonlinearElasticityandthePiolaTransformation 6
v 4. PreliminaryResultsforthePoisson-BoltzmannEquation 8
1
4.1. ThePoisson-BoltzmannequationwithPiolatransformation 8
7
3 4.2. RegularityandestimatesforthesingularsolutioncomponentG 9
1 4.3. Regularityandestimatesfortheregularlinearizedsolutioncomponentφr 12
. 4.4. Regularityandestimatesfortheregularnonlinearsolutioncomponentφr 15
1
0 5. AnElectrostaticForceModelandSomeEstimates 19
0 5.1. Thesurfaceforceduetochangingionicstrength 23
1 5.2. Thesurfaceforceduetoaddingalowdielectricconstantcavity 24
:
v 5.3. Thesurfaceforceduetoadditionalsingularcharges 26
i 5.4. Thesurfaceforceduetomolecularconformationalchange 27
X
5.5. Completeestimationoftheelectrostaticforces 30
r
a 6. MainResults: ExistenceofSolutionstotheCoupledSystem 30
7. VariationalPrincipleforExistenceand/orUniqueness 32
8. ConcludingRemarks 32
9. Acknowledgment 33
References 33
Date:July22,2007.
Key words and phrases. Macromolecular Conformational Change, Nonlinear Elasticity, Continuum
Modeling,Poisson-Boltzmannequation,ElectrostaticForce,CoupledSystem,FixedPoint.
MH was supported in part by NSF Awards 0411723, 022560 and 0511766, in part by DOE Awards
DE-FG02-04ER25620andDE-FG02-05ER25707,andinpartbyNIHAwardP41RR08605.
YZ and JAM were supported in part by the National Institutes of Health, the National Science Foun-
dation,theHowardHughesMedicalInstitute,theNationalBiomedicalComputingResource,theNational
ScienceFoundationCenterforTheoreticalBiologicalPhysics,theSanDiegoSupercomputingCenter,the
W.M.KeckFoundation,andAccelrys,Inc.
1
2 Y.C.ZHOU,M.HOLST,ANDJ.A.MCCAMMON
1. AN ELECTRO-ELASTIC MODEL OF CONFORMATIONAL CHANGE
A number of fundamental biological processes rely on the conformational change of
biomolecules and their assemblies. For instance, proteins may change their configu-
rations in order to undertake new functions, and molecules may not bind or optimally
bindtoeachothertoformnewfunctionalassemblieswithoutappropriateconformational
change at their interfaces or other spots away from binding sites. An understanding of
mechanisms involved in biomolecular conformational changes is therefore essential to
study structures, functions and their relations of macromolecules. Molecular dynamics
(MD)simulationshaveproventobeveryusefulinreproducingthedynamicsofatomistic
scalebytracingthetrajectoryofeachatominthesystem[37]. Despitetherapidprogress
made in the past decade mainly due the explosion of computer power and parallel com-
puting, it remains a significant challenge for MD to study large-scale conformational
changes occurring on time-scales beyond a microsecond [5]. Various coarse-grained
models and continuum mechanics models are developed in this perspective to comple-
ment the MD simulations and to provide computational tools that are not only able to
capture characteristics of the specific system, but also able to treat large length and time
scales. The prime coarse-grained approaches are the elastic network models, which de-
scribe the biomolecules to be beads, rods or domains connected by springs or hinges
accordingtothepre-analysisoftheirrigidityandtheconnectivity. Elasticnetworkmod-
els are usually combined with normal mode analysis (NMA) to extract the dominant
modes of motions, and these modes are then used to explore the structural dynamics at
reduced cost [9]. Continuum models do not depend on these rigidity or connectivity
analysis. Onthecontrary,therigidityofthestructureshallbeabletobederivedfromthe
results of the continuum simulations. Typical continuum models for biomolecular sim-
ulations include the elastic deformation of lipid bilayer membranes [35] and the gating
of mechanosensitive ion channels [34] induced by given external mechanical loads. It is
expected that with more comprehensive continuum models we will be able to simulate
the variation of the mechanical loads on the macromolecules with their conformational
change, and investigate the dynamics of molecules by coupling the loads and deforma-
tion. Thisarticletakesanimportantstepinthisdirectionbydescribingandanalyzingthe
first mathematical model for the interaction between the nonlinear elastic deformation
and the electrostatic potential field of macromolecules. Such coupled nonlinear models
have tremendous potential in the study of configuration changes and structural stability
oflargemacromoleculessuchasnucleicacids,ribosomesormicrotubulesduringvarious
electrostaticinteractions.
Ourmodelisdescribedbelow. LetΩ R3 beasmoothopendomainwhoseboundary
∈
is noted as ∂Ω; see Fig.(1). Let the space occupied by the flexible molecules Ω be a
mf
smoothsubdomainofΩ,whilethespaceoccupiedbytherigidmolecule(s)isdenotedby
Ω . Let the remaining space occupied by the aqueous solvent be Ω . The boundaries
mr s
of Ω and Ω are denoted by Γ and Γ , respectively. We assume that the distance
mf mr f r
betweenmolecularsurfacesand∂Ω
min x y : x Γ Γ ,y ∂Ω (1.1)
f r
{| − | ∈ ∪ ∈ }
is sufficiently large so that the Debye-Hu¨ckel approximation can be employed to de-
termine a highly accurate approximate boundary condition for the Poisson-Boltzmann
equation. TherearechargesatomslocatedinsideΩ andΩ ,andchangedmobileions
mf mr
in Ω . The electrostatic potential field generated by these charges induces electrostatic
s
forces on the molecules Ω and Ω . These forces will in turn cause the configuration
mf mr
ANONLINEARELASTICITYMODELOFCONFORMATIONALCHANGE 3
change of the molecules. We shall model this configuration rearrangement as an elas-
tic deformation in this study. Specifically, we will investigate the elastic deformation
of molecule Ω (which is originally in a free state and not subject to any net external
mf
force)inducedbyaddingmoleculeΩ andchangingmobilechargedensityinΩ . This
mr s
body deformation leads to the displacement of charges in Ω and the dielectric bound-
mf
aries, which simultaneously lead to change of the entire electrostatic potential field. It is
thereforeinterestingtoinvestigateifthedeformablemoleculeΩ hasafinalstablecon-
mf
figurationinresponsetotheappearanceofΩ andthechangeofmobilechargedensity.
mr
Ω
s Γ
! f
s Γr Ω!mf
m
Ω
mr
!
m
q
q i
j
∂Ω
FIGURE 1. Illustrationofmacromoleculesimmersedinaqueoussolventenvironment.
Withintheframeworkofanimplicitsolventmodelwhichtreatstheaqueoussolventin
Ω asastructure-lessdielectric,theelectrostaticpotentialfieldofthesystemisdescribed
s
bythePoisson-Boltzmannequation(PBE)
Nf+Nr
(cid:88)
((cid:15) φ)+κ2sinh(φ) = q δ(x )inΩ, (1.2)
i i
−∇· ∇
i
where δ(x ) is the Dirac distribution function at x , N + N is the number of singular
i i f r
charges of the system including the charges in Ω (i.e. N ) and Ω (i.e., N ). The di-
mf f mr r
electric constant (cid:15) and the modified Debye-Hu¨ckel parameter κ are piecewise constants
indomainsΩ ,Ω andΩ . Inparticular,κ = 0inΩ andΩ becauseitmodelsthe
mf mr s mf mr
free mobile ions which appear only in the solvent region Ω . The dielectric constant in
s
themoleculeandthatinthesolventaredenotedwith(cid:15) and(cid:15) ,respectively. Readersare
m s
referredto[26,27]fortheimportanceofthePoisson-Boltzmannequationinbiomolecu-
larelectrostaticinteractions,andto[2,3,4,28,29,30,31]forthemathematicalanalysis
aswellasvariousnumericalmethodsforthePoisson-Boltzmannequation.
The finite(large) deformation of molecules is essential to our coupled model, but can
not be described by a linear elasticity theory. We therefore describe the displacement
vectorfieldu(x)oftheflexiblemoleculeΩ withanonlinearelasticitymodel:
mf
div(T(u)) =f in Ω0
− b mf (1.3)
T(u) n =f on Γ0
· s f
wheref isthebodyforce,f isthesurfaceforceandT(u)isthesecondPiola-Kirchhoff
b s
stress tensor. In this study we assume the macromolecule is a continuum medium obey-
ingtheStVenant-Kirchhofflaw,andhenceitsstresstensorisgivenbythelinear(Hooke’s
law)stress-strainrelationforanisotropichomogeneousmedium:
T(u) = (I+ u)[λTr(E(u))I+2µE(u)].
∇
4 Y.C.ZHOU,M.HOLST,ANDJ.A.MCCAMMON
Hereλ > 0andµ > 0aretheLame´ constantsofthemedium,and
1
E(u) = ( uT + u+ uT u)
2 ∇ ∇ ∇ ∇
is the nonlinear strain tensor. The equation (1.3) is nonlinear due to the Piola transfor-
mation(I+ u)inT(u),andthequadraticterminthenonlinearstrainE(u). Thethird
∇
potential source of nonlinearity, namely a nonlinear stress-strain relation, is not consid-
eredhere;however,ourmethodsapplytothiscaseaswell.
It is noted that Eq.(1.3) is defined in the undeformed molecule body Ω0 with unde-
mf
formed boundary Γ0, while the Poisson-Boltzmann equation (1.2) holds true for real
f
deformed configurations. The deformed configuration is unknown before we solved
the coupled system. We therefore define a displacement-dependent mapping Φ(u)(x) :
Ω0 Ω and apply this mapping to the the Poisson-Boltzmann equation such that it can
→
also be analyzed on the undeformed molecular configuration. In Ω this map Φ(u)(x)
mf
is + u where is the identity mapping. A key technical tool in our work is that this
I I
mappingisthenharmonicallyextendedtoΩtoobtainthemaximumsmoothness. Apply
thismapping,thePoisson-Boltzmannequation1.2changestobe
Nf+Nr
(cid:88)
((cid:15)F(u) φ)+J(u)κ2sinh(φ) = J(u)q δ(Φ(x) Φ(x ))inΩ, (1.4)
i i
−∇· ∇ −
i
whereJ(u)istheJacobianofΦ(u)and
F(u) = ( Φ(u))−1J(u)( Φ(u))−T. (1.5)
∇ ∇
ThismatrixFiswelldefinedwheneverΦ(u)isaC1-diffeomorphism[6]. Thefunctions
in Eq.(1.4) should be interpreted as the compositions of respective functions in Eq.(1.2)
withmappingΦ(x),i.e.,φ(x) = φ(Φ(x)),(cid:15)(x) = (cid:15)(Φ(x))andκ(x) = κ(Φ(x)).
In this paper, we shall analyze the existence of the coupled solution of the elastic-
ity equation (1.3) and the transformed Poisson-Boltzmann equation (1.4). These two
equations are coupled through displacement mapping Φ(u) in the Poisson-Boltzmann
equation and the electrostatic forces to be defined later. The solution of this coupled
system represents the equilibrium between the elastic stress of the biomolecule and the
electrostatic forces to which the biomolecule is subjected. The existence, the unique-
ness and the W2,p-regularity of the elasticity solution have already been established by
Grandmont [6] in studying the coupling of elastic deformation and the Navier-Stokes
equations; thus in this work we shall focus on the solution to the transformed Poisson-
Boltzmann equation and to the coupled system. We shall define a mapping S from an
appropriate space X of displacement field u into itself, and seek the fixed-point of this
p
map. This fixed-point, if it exists, will be the solution of the coupled system. A critical
step in defining S is the harmonic extension of the Piola transformation from Ω to
mf
Ω and R3. The regularity of the Piola transformation determines not only the existence
of the solution to the transformed Poisson-Boltzmann equation, but also the existence
of the solution to the coupled system. Because most of our analysis will be carried out
ontheundeformedconfigurationwewillstilluseΩ ,Ω ,Ω ,Γ ,Γ todenotetheun-
mf mr s f r
deformed configurations of molecules, the solvent and the molecular interfaces, unless
otherwisespecified.
The paper is organized as follows. In Section 2 we review a fundamental result
concerning the piecewise W2,p-regularity of the solutions to elliptic equations in non-
divergence form and with discontinuous coefficients. The nonlinear elasticity equation
will be discussed in Section 3, where the major results from [6] are presented without
ANONLINEARELASTICITYMODELOFCONFORMATIONALCHANGE 5
proof. The Piola transformation will be defined, harmonically extended, and then an-
alyzed. In Section 4 we will prove the existence and uniqueness of the solution to the
Piola-transformed Poisson-Boltzmann equation, generalizing the results in [2] for the
un-transformed case. Both L∞ and W2,p estimates will be given for the electrostatic
potential in the solvent region, again generalizing results in [2]. We will then define the
electrostaticforcesandestimatetheseforcesbydecomposingthemintocomponentscor-
responding to four independent perturbation steps. The estimates of these components
areobtainedseparatelyandthefinalestimateofthesurfaceforceisassembledfromthese
individual estimates. The coupled system will be finally considered in Section 6 where
the mapping S will be defined, and the main result of the paper will be established by
applying a fixed-point theorem on this map to give the existence of a solution of the
coupledsystem.
2. NOTATION AND SOME BASIC ESTIMATES
In what follows Wk,p( ) will denote the standard Sobolev space on an open domain
D
, where can be Ω,Ω or Ω . While solutions of the Poisson-Boltzmann have low
m s
D D
global regularity in Ω, we will need to explore and exploit the optimal regularity of the
solution in any subdomain of Ω. For this purpose, we define 2,p(Ω) = W2,p(Ω ) (cid:117)
m
W2,p(Ω ) where (cid:117) is the direct sum. Every function φ 2,pWcan be written as φ(x) =
s
∈ W
φ (x)+φ (x)whereφ (x) W2,p(Ω ),φ (x) W2,p(Ω ),andhasanorm
m s m m s s
∈ ∈
φ = φ + φ . (2.1)
(cid:107) (cid:107)W2,p (cid:107) m(cid:107)W2,p(Ωm) (cid:107) s(cid:107)W2,p(Ωs)
Similarly, we define a class of functions = (Ω) which are continuous in either sub-
C C
domain and may have finite jump on the interface, i.e., a function a is given by
∈ C
a = a +a where a C(Ω ),a C(Ω ) are continuous functions in their respec-
m s m m s s
∈ ∈
tivedomains. Thenormin isdefinedby
C
a = a + a .
(cid:107) (cid:107)C (cid:107) m(cid:107)C(Ωm) (cid:107) s(cid:107)C(Ωs)
We recall two important results. The first is a technical lemma which will be used for
theestimationoftheproductoftwoW1,p functions;thisissometimescalledtheBanach
algebraproperty.
Lemma 2.1. Let 3 < p < ,1 q p be two real numbers. Let Ω be a domain in R3.
∞ ≤ ≤
Let u W1,p(Ω),v W1,q(Ω), then their product uv belongs to W1,q, and there exists
∈ ∈
aconstantC suchthat
uv C u v .
W1,q(Ω) W1,p(Ω) W1,q(Ω)
(cid:107) (cid:107) ≤ (cid:107) (cid:107) (cid:107) (cid:107)
For the proof of this lemma we refer to [1]. In this paper we will apply Lemma (2.1)
to the case with p = q. The second result is a theorem concerning the Lp estimate of
ellipticequationswithdiscontinuouscoefficients.
Theorem 2.2. Let Ω and Ω Ω be bounded domains of R3 with smooth boundaries
1
⊂⊂
∂ΩandΓ. LetΩ = (Ω Γ)andΩ = Ω Ω . LetAbeasecondorderellipticoperator
1 1 2 1
∪ \
suchthat
(cid:40)
(A u)(x) x Ω (cid:88)
(Au)(x) = 1 ∈ 1 , whereA = a (x)Dk.
i ik
(A u)(x) x Ω
2 2
∈ k≤2
6 Y.C.ZHOU,M.HOLST,ANDJ.A.MCCAMMON
Thenthereexistsauniquesolutionu 2,p fortheinterfaceproblem
∈ W
(Au)(x) = f inΩ
[u] = u u = 0onΓ
2 1
−
[Bu ] = B u n B u n = honΓ
n 2 2 1 1
∇ · − ∇ ·
u = g on∂Ω
providingthata (Ω),B C(Γ),f Lp(Ω),g W2−1/p,p(∂Ω),h W1−1/p,p(Γ),
ik i
∈ C ∈ ∈ ∈ ∈
wherenistheoutsidenormaltoΩ . Moreover,thefollowingestimateholdstrue
1
(cid:0) (cid:1)
u K f + h + g + u , (2.2)
W2,p(Ω) Lp(Ω) W1−1/p,p(Γ) W2−1/p,p(∂Ω) Lp(Ω)
(cid:107) (cid:107) ≤ (cid:107) (cid:107) (cid:107) (cid:107) (cid:107) (cid:107) (cid:107) (cid:107)
wheretheconstantK dependsonlyonΩ,Ω ,Ω ,pandthemodulusofcontinuityofA.
1 2
Theorem (2.2) is fundamental to various results about elliptic equations with discon-
tinuouscoefficients;Forexample,theglobalH1 regularityandH2 estimatesofBabu¨ska
[20], the finite element approximation of Chen et al. [21], a prior estimates for second-
orderellipticinterfaceproblems[22],thesolutiontheoryandestimatesforthenonlinear
Poisson-Boltzmann equation [2], and the continuous and discrete a priori L∞ estimates
for the Poisson-Boltzmann equation along with a quasi-optimal a priori error estimate
for Galerkin methods [2] applied to the Poisson-Boltzmann equation. For the proof of
Theorem (2.2) and the more general conclusions for high-order elliptic equations with
high-orderinterfaceconditionswereferto[23,24].
3. NONLINEAR ELASTICITY AND THE PIOLA TRANSFORMATION
We first state a theorem concerning the existence, uniqueness, regularity and the esti-
mationofthesolutiontothenonlinearelasticityequation[6]:
Theorem3.1. Letthebodyforcef Lp(Ω )andthesurfaceforcef W1−1/p,p(Γ ),
b mf s f
∈ ∈
where 3 < p < . There exists a neighborhood of 0 in Lp(Ω ) W1−1/p,p(Γ )
mf f
∞ ×
such that if (f ,f ) belongs to this neighborhood then there exists a unique solution u
b s
∈
W2,p(Ω ) W1,p (Ω )of
mf ∩ 0,Γf0 mf
div(T(u)) =f in Ω ,
b mf
−
T(u)n =f on Γ Γ ,
s f f0
\
u =0on Γ , (3.1)
f0
(cid:90)
(I+ u)J(u)(I+ u)−T n =3 Ω ,
mf
∇ ∇ · | |
Γ
f
where Γ is a subset of Γ equipped with homogeneous Dirichlet boundary condition,
f0 f
I is the unit matrix. The last equation represents the incompressibility condition of the
elastic deformation. Moreover, the solution can be estimated with respect to the force
data:
u C( f + f ). (3.2)
(cid:107) (cid:107)W2,p(Ωmf) ≤ (cid:107) b(cid:107)Lp(Ωmf) (cid:107) s(cid:107)W1−1/p,p(Γf)
Proof. See[7]. (cid:3) (cid:3)
Remark3.2. Itisnoticedthatu C1,1−3/p(Ω )becauseofthecontinuousembedding
mf
∈
ofW2,p(Ω )inC1,1−3/p(Ω )forp > 3.
mf mf
Thedisplacementfieldu(x)solvedfromEq.(3.1)naturallydefinesamappingΦ(u) =
+uinΩ where istheidentitymapping. ThismappingΦ(u)(x)hastobeappropri-
mf
I I
atelyextendedintoR3 Ω toyieldaglobaltransformationforthePoisson-Boltzmann
mf
\
ANONLINEARELASTICITYMODELOFCONFORMATIONALCHANGE 7
equation. Itiscriticalinwhatfollowsthatthisextensionhasvariousfavorableproperties,
whichleadsustodefineaglobalmappingbyharmonicextension:
(cid:40)
+u x Ω
mf
Φ(u) = I ∈ (3.3)
+w otherwise
I
wherew solves
∆w =0inR3 Ω ,
mf
\ (3.4)
w =uonΓ .
f
The following crucial lemma concerns the regularity of Φ(u) and the invertibility of
Φ(u):
∇
Lemma3.3. LetΦ(u)bedefinedinEq.(3.3),wehave
(a) Φ(u) W2,p(Ω )andΦ(u) C∞(R3 Ω ).
mf mf
∈ ∈ \
(b) There exists a constant M > 0 such that for all u M, Φ(u) is an
W2,p(Ω)
(cid:107) (cid:107) ≤ ∇
invertiblematrixinW1,p(Ω )andinC∞(R3 Ω ).
mf mf
(c) Undercondition(b),Φ(x)isone-to-oneonR3,\andisaC1-diffeomorphismfrom
Ω to Φ(u)(Ω ), and is a C∞-diffeomorphism from R3 Ω to Φ(u)(R3
mf mf mf
\ \
Ω ).
mf
Proof. That Φ(u) W2,p(Ω ) follows directly from its definition. Also, of Φ(u)
mf
∈ ∈
C∞(R3 Ω ) since Φ(u) = + w while w is harmonic hence analytical in Φ(u)
mf
\ I ∈
C∞(R3 Ω )becauseitisthesolutionoftheLaplaceequation(3.4). Fortheinvertibility
mf
\
of Φ(u) in W1,p(Ω ) we refer to Lemma 2 in [6] or Theorem 5.5.1 in [7], which says
mf
thatifau Ω isdifferentiableand
mf
∈
u(x) < C
|∇ |
for some constant depending on Ω , then Φ(u) = I + u > 0 x Ω and
mf mf
∇ ∇ ∀ ∈
I+ u is injective on Ω . The invertibility of Φ(u) therefore follows from the facts
mf
∇ ∇
thatu C1,1−3/p(Ω )suchthatforsufficientlysmallM
mf
∈
u u C u = C M C.
|∇ | ≤ (cid:107) (cid:107)C1,1−3/p(Ωmf) ≤ 1(cid:107) (cid:107)W2,p(Ωmf) 1 ≤
To prove the invertibility of Φ(u) = I + w in R3 Ω we notice the following
mf
∇ ∇ \
estimateforthefirstderivativeofthesolutiontoLaplaceequation[19]:
w w C u C M C,
|∇ | ≤ (cid:107) (cid:107)C1(R3\Ωmf) ≤ 2(cid:107) (cid:107)C1,1−3/p(Γf) ≤ 2 ≤
ThereforeifM ischosensuchthat
C
M (3.5)
≤ max C1,C2
{ }
Φ(u)isaninvertiblematrixinR3. (cid:3) (cid:3)
∇
Remark 3.4. It follows from Lemma (3.3) that the matrix F(u) in Eq.(1.5) is well-
defined,symmetricandpositivedefinite. Moreprecisely,wehavethatthemapsF(u)(x)
∈
C0,1−3/p(Ω ) and F(u)(x) C∞(R3 Ω ). On the other hand, as a mapping from
mf mf
∈ \
u W2,p(Ω ) to F(u) C∞(R3 Ω ), F(u) is infinitely differentiable with respect
mf mf
∈ ∈ \
to u. In all what follows we will write F(u) and J(u) as F and J only, keeping in mind
thattheyareudependent.
8 Y.C.ZHOU,M.HOLST,ANDJ.A.MCCAMMON
4. PRELIMINARY RESULTS FOR THE POISSON-BOLTZMANN EQUATION
4.1. ThePoisson-BoltzmannequationwithPiolatransformation. Therigorousanal-
ysis and numerical approximation of solutions to the Poisson-Boltzmann equation (1.2)
or its transformed version (1.4) are generally subject to three major difficulties: 1) the
singular charge distribution, 2) the discontinuous dielectric constant on the molecular
surface and 3) the strong exponential nonlinearities. However, it was recently demon-
strated [2] that as far as the untransformed Poisson-Boltzmann equation (1.2) is con-
cerned, some of these difficulties can be side-stepped by individually considering the
singular and the regular components of the solution. Specifically, the potential solution
isdecomposedtobe
φ = G+φr = G+φl +φn (4.1)
wherethesingularcomponent
(cid:88) q
i
G =
(cid:15) x x
m i
i | − |
isthesolutionofthePoissonequation
N
(cid:88)
((cid:15) G) = ρ := q δ(x ) in R3; (4.2)
m f i i
−∇· ∇
i
whileφl isthelinearcomponentoftheelectrostaticpotentialwhichsatisfies
((cid:15) φl) = (((cid:15) (cid:15) ) G) in Ω,
m
−∇· ∇ −∇· − ∇ (4.3)
φl =g G on ∂Ω,
−
andthenonlinearcomponentφn solves
((cid:15) φn)+κ2sinh(φn +φl +G) =0 in Ω,
−∇· ∇ (4.4)
φn =0 on ∂Ω,
where
(cid:88)N e−κ|x−xi|
g = q (4.5)
i
(cid:15) x x
s i
i=1 | − |
is the boundary condition of the complete Poisson-Boltzmann equation (1.2). Such a
decomposition scheme removes the point charge singularity from the original Poisson-
Boltzmannanditwasshownin[2]thattheregularcomponentoftheelectrostaticpoten-
tial φr = φl +φn belongs to H1(Ω) although the entire solution G+φr does not. The
most prominent advantage of this decomposition lies in the fact that the regular compo-
nent represents the reaction potential field of the system, which can be directly used to
compute the solvation energy and other associated important properties of the system.
It is not necessary to solve the Poisson-Boltzmann equation twice, once with uniform
vacuum dielectric constant and vanishing ionic strength and the other with real physical
conditions,toobtainthereactionfield[29]. Astobeshownlateron,theidentificationof
this regular potential component as the reaction field also facilitates the analysis and the
computationoftheelectrostaticforces.
Applying the similar decomposition to the transformed Poisson-Boltzmann equation
wegetanequationforthesingularcomponentG:
((cid:15) F G) = Jρ in R3, (4.6)
m f
−∇· ∇
ANONLINEARELASTICITYMODELOFCONFORMATIONALCHANGE 9
andanequationfortheregularcomponentφr:
((cid:15)F φr)+Jκ2sinh(φr +G) = (((cid:15) (cid:15) )F G) in Ω,
m
−∇· ∇ ∇· − ∇ (4.7)
φr =g G on ∂Ω.
−
We shall prove the existence of φr in Eq.(4.7) and give its L∞ bounds by individually
consideringtheequationforthelinearcomponentφl:
((cid:15)F φl) = (((cid:15) (cid:15) )F G) in Ω,
m
−∇· ∇ ∇· − ∇ (4.8)
φl =g G on ∂Ω,
−
andtheequationforthenonlinearcomponentφn:
((cid:15)F φn)+Jκ2sinh(φn +φl +G) =0 in Ω,
−∇· ∇ (4.9)
φn =0 on ∂Ω.
As mentioned above, the functions G,φl,φn,ρf and κ in Eqs.(4.6) through (4.9) shall
be interpreted as the compositions of the corresponding entries of these functions in
untransformedequations(4.2)through(4.4)withthePiolatransformationΦ(x),i.e.,g =
g(Φ(x)),G = G(Φ(x)),φl = φl(Φ(x)),φn = φn(Φ(x)),ρf = ρf(Φ(x)),κ = κ(Φ(x)).
4.2. RegularityandestimatesforthesingularsolutioncomponentG. Wefirststudy
the Eq. (4.6) for the singular component of electrostatic potential. We remark that the
linearandnonlinearPBequationshavethesamesingularcomponentoftheelectrostatic
potential. ThesolutionofthissingularcomponentistheGreen’sfunctionfortheelliptic
operatorLdefinedby
Lu = ((cid:15) F u). (4.10)
m
−∇· ∇
Weshallusethefollowingtheorem[11]concerningtheregularityandtheestimateofthe
Green’sfunction:
Theorem4.1. LetΩbeanopensetinR3. Supposetheellipticoperator
n
(cid:88) ∂ ∂u
Lu = (a )
ij
∂x ∂x
j i
i,j=1
isuniformlyellipticandbounded,whilethecoefficientsa satisfying
ij
a (x) a (y) ω( x y )
ij ij
| − | ≤ | − |
foranyx,y Ω,andthenon-decreasingfunctionω(x)satisfies
∈
ω(2t) Kω(t)forsomeK > 0andallt > 0,
≤
(cid:90)
ω(t)
dt < .
t ∞
R
ThenforthecorrespondingGreen’sfunctionGthefollowingsixinequalitiesaretruefor
anyx,y Ω:
∈
(a) G(x,y) K x y −1,
≤ | − |
(b) G(x,y) Kδ(x) x y −2.
≤ | − |
(b) G(x,y) Kδ(x)δ(y) x y −3.
≤ | − |
(d) G(x,y) K x y −2.
x
|∇ | ≤ | − |
(e) G(x,y) Kδ(y) x y −3.
y
|∇ | ≤ | − |
(f) G(x,y) K x y −3.
x y
|∇ ∇ | ≤ | − |
whereδ(y) = dist(y,∂Ω)andthegeneralconstantK = K(a ,ω,Ω).
ij
10 Y.C.ZHOU,M.HOLST,ANDJ.A.MCCAMMON
FromthistheoremwecanderivetheregularityoftheGreen’sfunctionoftheoperator
(4.10). Indeed, by Sobolev embedding (cid:15) F C0,1−3/p(R3), therefore it satisfies the
m
∈
conditions on a in this theorem provided that ω(t) = Kt3/p. We then conclude that the
ij
singular component of the electrostatic potential G W1,∞(Ω B (x )). On the other
r i
∈ \
hand, from Eq. (4.6) we know that G(Φ(u)(x))/J(x ) itself is the Green’s function of
i
operator (4.10) if F is generated by the Piola transformation according to (1.5) and J
is the corresponding Jacobian. Thus the Green’s function of differential operator (4.10)
belongstoW2,p(Ω B (x ))sinceitisthecompositionoftheGreen’sfunctionofLaplace
r i
\
operator,whichisofC∞(Ω B (x )),andthePiolatransformation,whichisofW2,p(Ω).
r i
\
Higher regularity of G in Ω can be derived thanks to the harmonic extension of u to
s
R3 Ω . In particular, because all charges are located in Ω and Ω the Poisson
mf mf mr
\
equation(4.6)appearsaLaplaceequation
((cid:15)F G) = 0 in Ω ,
s
∇ ∇
henceG(x) C∞(Ω ),sinceΩ isasmoothopendomainandF C∞(Ω ).
s s s
∈ ∈
InadditiontotheregularityoftheGreen’sfunction,wehavefollowingestimatesofG
withrespecttoFandJ.
Lemma4.2. ForanygivenmoleculetheGreen’sfunctionGofoperator(4.10)hasesti-
mates
(a) G C J .
(cid:107) (cid:107)L∞(Ωs) ≤ (cid:107) (cid:107)L∞(Ω)
(b) G C J .
(cid:107)∇ (cid:107)L∞(Ωs) ≤ (cid:107) (cid:107)L∞(Ω)
If in addtion F I C , J 1 C for some constant C and C ,
W1,p(Ω) f W1,p(Ω) J f J
(cid:107) − (cid:107) ≤ (cid:107) − (cid:107) ≤
then
(c) G C G .
(cid:107) (cid:107)Lp(∂Ω) ≤ (cid:107) (cid:107)L∞(Ωs)
(d) g Φ C g .
(cid:107) ◦ (cid:107)W2−1/p,p(∂Ω) ≤ g(cid:107) (cid:107)W2,p(Ωs)
(e) g Φ G C g +C G .
(cid:107) ◦ − (cid:107)W2−1/p,p(∂Ω) ≤ g(cid:107) (cid:107)W2,p(Ωs) G(cid:107) (cid:107)L∞(Ωs)
(f) F G C G forsomesetΩ(cid:48).
(cid:107) ∇ (cid:107)W1−1/p,p(Γ) ≤ Γ(cid:107) (cid:107)L∞(Ω(cid:48)s) s
Proof. This J is well defined since J is uniformly continuous in Ω . To prove
(cid:107) (cid:107)L∞(Ωs) s
(a)and(b)wedefineq = max q and
max i
{| |}
K K
G (x,x ) = , G (x,x ) =
(cid:107)∇x i i (cid:107)L∞(Ωs) δ2 (cid:107) i i (cid:107)L∞(Ωs) δ
whereδisthesmallestdistancebetweenx ∂Ωandsingularchargesatx . Thissmallest
i
∈
distance is related to the radii of atoms used in defining the molecular surface. In the
sense of Connolly’s molecular surface, δ is simply the smallest van der Waals radius of
the atoms which have contact surface [33]. We can therefore bound G and its gradient
with
(cid:88)
G = Jq G Nq J G
(cid:107) (cid:107)L∞(Ωs) (cid:107) i i(cid:107)L∞(Ωs) ≤ max(cid:107) (cid:107)L∞(Ωs)(cid:107) i(cid:107)L∞(Ωs)
i
J NKq
L∞(Ω) max
= (cid:107) (cid:107) , (4.11)
δ
(cid:88)
G = Jq G Nq J G
(cid:107)∇ (cid:107)L∞(Ωs) (cid:107) i∇ i(cid:107)L∞(Ωs) ≤ max(cid:107) (cid:107)L∞(Ωs)(cid:107)∇x i(cid:107)L∞(Ωs)
i
J NKq
L∞(Ω) max
= (cid:107) (cid:107) , (4.12)
δ2
whereN isthetotalnumberofsingularchargesand J isthemaximumJacobian
(cid:107) (cid:107)L∞(Ωs)
onΓ.