Schlick T. Molecular Modeling and Simulation: An Interdisciplinary Guide

Подождите немного. Документ загружается.

330 10. Nonbonded Computations

are determined by an interaction list). These multipole expansions for the clusters

are then used to compute interactions between distant clusters of particles, with

nearby interactions computed directly.

This clustering is at the heart of such hierarchical methods: the computation for

the interactions between particles is recast as work between clusters containing

groups of particles.

O(N log N) Work

The recursive computation described above is completed in about log

N steps

(levels), leading to total work of O(N log N).However,theconstant associated

with (N log N) depends on the number of operations required to form the expan-

sions. This cost of the local expansions depends on the number of coefﬁcients p

used (like the number of Fourier terms in FFTs). This number can be speciﬁed

based on the desired accuracy, 

acc

, for the resulting potential.

It can be shown that p can be related to 

acc

p =log

√

(1/

acc

)

at sufﬁciently large separations. The corresponding total work is approximately

200 Np

log

N [481].

Fast Multipole Machinery (O(N))

Unfortunately, these estimates do not produce substantial speedups in practice

with respect to direct calculations (considering also the additional bookkeeping

work required for multipole schemes), since p ≈ 20 for 7 digits of accuracy.

Signiﬁcant speedups are achieved only for very large system sizes.

Independently of the Barnes and Hut scheme, the fast multipole method was

also developed in the mid 1980s by Rokhlin [1065]. Greengard and Rokhlin then

made a seminal contribution [479] by developing further mathematical machin-

ery that exploits the smoothness of the far-ﬁeld potential and works with the

tree codes. Relying on translation operators that act on the distant (multipole)

as well as local expansions (harmonics), they lumped and converted expansions

associated with several clusters into local expansions. The cost of translating the

multipole-to-local expansions is reduced via fast schemes for application of ro-

tation matrices. This combined clever mathematical and numerical machinery

yields an asymptotically O(N) method. The linear constant depends sensitively

on the practical algorithmic implementation.

10.5.3 Expansion in Spherical Coordinates

The multipole expansion is most conveniently written in spherical coordinates.

A point x in 3D space can be represented by the triplet {r, θ, φ} instead of the

Cartesian triplet {x, y, z},where:

r =

+ y

+ z

,θ=cos

−1

(z/r),φ=tan

−1

(y/z).

10.5. The Multipole Method 331

Consider N charges located at points represented in spherical coordinates as

{ρ

,α

,β

} lying inside a sphere of radius a, i.e., with |ρ

| <afor all i.Then

for any 3D point x ≡{r, θ, φ} outside that sphere, i.e., with r>a, the potential

at x can be written in terms of spherical harmonics functions, {Y

}, solutions

of the Laplace equation in spherical coordinates, as follows:

Φ(x)=

∞



n=0



m=−n

n+1

(θ, φ) . (10.49)

The functions Y

(θ, φ) are called the spherical harmonics of degree −(n +1),

or multipoles, and the coefﬁcients M

are known as moments of the expansion:



i=1

−m

(α

,β

) . (10.50)

The spherical harmonics functions can be expressed in terms of partial

derivatives of 1/r from the following relations [477]:

(θ, φ)=

⎧

⎪

⎨

⎪

⎩

n+1

∂

∂z





m =0

n+1



∂

∂x

+ i

∂

∂y





∂

∂z



(n−m)





m =1, 2, ...

n+1

−m



∂

∂x

− i

∂

∂y



−m



∂

∂z



(n+m)





m = −1, −2, ...

(10.51)

where

(−1)

(n − m)!(n + m)!

. (10.52)

The spherical harmonic functions {Y

(θ, φ)} are also related to the associated

Legendre polynomials of degree n, {P

},deﬁnedas

(x)=(−1)

(1 −x

)

m/2

(x)

(θ, φ) ≡

(n −|m|)!

(n + |m|)!

|m|

(cos θ)exp(imφ) .

The expansion power p is chosen to achieve the desired accuracy so that the

remainder



∞



n=p+1



m=−n

n+1

(θ, φ)



≤



i=1

r − a





p+1

= 

acc

. (10.53)

332 10. Nonbonded Computations

Thus to achieve greater accuracy, we can increase the expansion power p

or decrease the radius a. For example, if r =2a, the remainder above is

−(p+1)

/a)



|. Setting p =log

(

−1

acc

) yields an accuracy of 

acc

relative



|/a.

The spherical harmonics functions are believed to provide a more efﬁcient basis

than Cartesian expansions, since they involve p

coefﬁcients rather than p

10.5.4 Biomolecular Implementations

Thorough comparisons of the balance among speed, accuracy, and scalability

of Ewald and fast multipole schemes are ongoing areas of research. The fast

multipole approach was ﬁrst implemented for biomolecular dynamics in the

early 1990s [142, 1176], and later parallelized to achieve good speedup on a

large number of processors [140, 141, 1133], both on workstation clusters and

on supercomputers like Cray T3D/T3E. A periodic version of the fast multipole

method relies on machinery similar to that used for the periodic Ewald summa-

tion [1136], namely obtaining a convergent sum by masking the original sum by

Gaussians.

Implementations in 3D of the periodic fast multipole method were adapted for,

and applied to, molecular dynamics by Board and co-workers in the late 1990s

[693, 1133] and by Figueirido et al. [399], as well as compared to the PME

alternative.

In 3D applications, implementors have found it particularly difﬁcult to achieve

good performance with high accuracy. For example, in [399] it has been reported

that the multipole approach is three times slower than PME for about 6000 atoms

and slightly slower than PME for about 22,000 atoms; it is suggested, however,

that fast multipoles will be competitive for larger sizes.

In serial test implementations, Board and collaborators found that the periodic

version of the fast multipole method is still signiﬁcantly slower than PME by

roughly a factor of two for a large water system of 71,496 atoms [1133]; this

performance involves modest accuracy (4–5 digits of accuracy in the potential

and 3–4 in the corresponding force) and required 8 multipole terms. This per-

formance result may be explained by the very fast, ‘custom’ version of PME of

Board and co-workers which uses a table lookup for the erfc function and yields

equivalent run times for calculations involving spherical cutoffs of 10–12

A. Work

for distributed PME versus distributed parallel multipole tree code, run at the

modest accuracy, still shows PME to be faster on a small number of processors

(Figure 10.8).

The newer version of the fast multipole method extends the 2D scheme of trans-

lation operators acting on harmonic functions to produce a new diagonal form for

these operators [481]. This approach produces a more modest break-even sys-

tem size in 3D when compared to direct calculations, such as 2000 particles for

single-precision accuracy or 5000 for 10-digit accuracy. Application of this fast

multipole method to biomolecular dynamics should be forthcoming.

10.6. Continuum Solvation 333

10.5.5 Other Variants

Besides these commonly-used fast multipole and Ewald methods, scientists have

also developed variants as well as new methods. Summations based on multi-

grid techniques developed for partial differential equations, as implemented for

the fast evaluation of integral transforms [164], were extended to molecular elec-

trostatic potentials [1088,1198]. A variable-order extension of Appel’s algorithm

in Cartesian coordinates was also reported, with adaptive tree structure and error

control [342,343].

10.6 Continuum Solvation

10.6.1 Need for Simpliﬁcation!

The fast electrostatic techniques described above have been designed and are nec-

essary for ‘deluxe’ models of macromolecular systems — very detailed solvated

biomolecular representations involving thousands of explicit water molecules sur-

rounding the biomolecule or biomolecules. This large size mimics the cellular

environment and is needed to reproduce the bulk properties of water as well as

the effects of local solvation on important associations such as those between two

proteins, proteins and ligands, and proteins and nucleic acids.

While such elaborate models can yield invaluable information on the rich com-

plexity of biomolecular environments, they are clearly computationally-intensive

and, unfortunately, free neither from approximations nor artifacts. As described

previously, artiﬁcial periodicity can lead to noncanceling errors associated with

the Coulomb and solvation contributions to the electrostatic free energy [580](see

also below). Trajectories are also highly sensitive to various force-ﬁeld approxi-

mations, propagation protocols, and modeling choices (e.g., positions of ions and

water molecules); see discussion in the chapters on molecular dynamics.

Implicit solvent approaches reduce the number of degrees of freedom dras-

tically and account for them in an average sense, in the form of solvation free

energy estimates. Such approaches can be effective, especially when combined

with coarse grained models of molecular systems. Practical methods that involve

clever engineering constructs of various approximations have been developed by

many researchers, including McCammon, Case, Karplus, Roux, Honig, Truh-

lar, and many others, as recently reviewed [223, 270, 377, 587]. However, Chen

and Brooks caution that current surface-area based nonpolar models have signif-

icant limitations and thus could beneﬁt by the incorporation of several nonpolar

solvation aspects [223].

Thorough introductions to implicit solvation models can be found in [758,

1073], the latter with emphasis on quantum-mechanical based approaches with a

literature survey of relevant works. See [223] for a more recent review including

many technical details, Cramer’s textbook [270] for a thorough treatment, and

[376,1190,1263] for reviews of Poisson-Boltzmann and generalized Born models.

334 10. Nonbonded Computations

10.6.2 Potential of Mean Force

Implicit solvent models based on continuum electrostatics treatments form an

alternative to this deluxe approach to modeling [813, 1073, 1190, 1263]. These

methods essentially approximate a potential of mean force [648] for the solvated

biomolecular environment and can provide a broad, qualitative and quantitative

view of overall biomolecular structure, recognition, and functional properties.

Balancing Biophysics with Numerics

The potential of mean force is an effective free energy potential that depends

on state variables like temperature and pressure [648]. It augments the potential

energy used to describe the intramolecular interactions by indirectly incorporat-

ing the effects of rearrangement in the medium on the biomolecule (solute). This

effect is important to consider since water molecules and ions in the solute en-

vironment perturb, in a cooperative manner, the pairwise forces acting on pairs

of solute particles through their own locations. The basic idea is to construct a

function that captures in an average sense solute conﬁgurations without explicit

consideration of the solvent degrees of freedom. Additionally, hydrodynamicscan

be applied to provide continuum approximations to the effects on the solute of

the motions of the solvent molecules. In this way, the model complexity can be

compromised with physical reproducibility/reliability.

Electrostatic and Non-Electrostatic Components

When such potentials of mean force are constructed, they are generally sep-

arated into two components: electrostatic and non-electrostatic (or nonpolar)

interactions [580, 1073]. This decomposition is useful for dissecting the differ-

ent microscopic factors that inﬂuence molecular conformations and solvation.

The electrostatic component can also be related to various continuum electrostatic

approximations (see below).

Variations

There are many ways to construct the effective potential of mean force for this

broader structural and functional description, both empirically and theoretically.

Empirical constructs lump all effects into an ‘information based’ or ‘statistical’

potential, derived, for example, from observed protein structures in solution [912].

Other approaches involve approximations of integral equations [1014], free

energy simulations, solvent-accessible surface models, various combinations

of implicit/explicit solvent representations [150, 785, 1314], generalized Born

models [96, 1273], solutions based on classical continuum electrostatics (i.e.,

Poisson-Boltzmann equation) [565] and quantum mechanics [758], and phe-

nomenological dynamic models such as Langevin and Brownian dynamics (BD)

[813]. Generalized Born models, in particular, construct a dielectric cavity sur-

rounding the molecular charges to approximate solvation and salt-screening

effects. In all cases, the dielectric constant must be chosen with care [1146].

10.6. Continuum Solvation 335

A good review on implicit solvent models can be found in [223]. See also

the special volume in Biophysical Chemistry, volume 78 (5 April 1999), guest

edited by Benoˆıt Roux and Thomas Simonson, with updates in [1190,1263]. The

study in [1439] shows that implicit solvent models yield reasonable agreement in

terms of protein/ligand binding free energies when compared to their much more

computational-costly explicit-solvent models. Applications to membrane systems

[1263] reveal the progress of modeling such heterogeneous materials by implicit

solvation models but also some computational problems (e.g., artiﬁcial period-

icity imposed by Ewald sums and related methods) that require further work to

characterize and resolve.

Below, we sketch two approaches, based on Langevin dynamics and continuum

electrostatics.

10.6.3 Stochastic Dynamics

The theory of hydrodynamics has a long history in molecular modeling, orig-

inating from the study of liquids, polymers, and simple molecular reactions

[184, 966, 1072, 1407, 1456]. These theories have been applied to macromolec-

ular dynamics via Langevin and Brownian dynamics simulations that generate

stochastic trajectories, so called because the governing dynamic equations in-

clude stochastic forces that mimic solvent effects, in addition to the systematic

force (negative gradient of the potential energy).

See classic statistical mechanics texts such as [853] for extensive background,

and Chapter 14 of this text for a brief discussion of generating stochastic trajec-

tories by Langevin and Brownian dynamics algorithms; a thesis by Hongmei Jian

[607] nicely summarizes the theory and numerical applications of Langevin and

Brownian dynamics to long DNA molecules.

In the stochastic treatment, the inﬂuence of solvent particles on the solute is

incorporated through additional frictional and random terms in a manner consis-

tent with physical laws regarding equilibrium and nonequilibrium processes (e.g.,

equilibrium conformation distributions, ﬂuctuation/dissipation theorem) [853].

Applications of Langevin and Brownian dynamics simulations have been partic-

ularly successful for macroscopic models of biomolecules, such as long DNA of

thousands of base pairs [107, 233, 608, 1106]. In such applications, polymer the-

ory has been used to guide parameterization (for hydrodynamic radii, timesteps,

etc.) through governing macroscopic polymer properties [410] such as persistence

length, radius of gyration, and diffusion constants. Protein applications include

long-timescale enzyme catalysis events such as the loop opening/closing motion

in triosephosphate isomerase (TIM) [298,299,1321,1322].

The Langevin Equation

In the simplest form of the Langevin equation, the friction kernel is taken to be

space and time independent for each particle, and the inﬂuence of the environment

on the systematic, internal force is represented in an average sense. Thus, explicit

336 10. Nonbonded Computations

hydrodynamic interactions are ignored, and the internal force is augmented by a

frictional term, proportional to the velocity, and a random force R, which crudely

mimics molecular collisions and viscosity in the realistic cellular environment:

= −∇E(x) − Mγ

+ R(t). (10.54)

Here E is the potential energy governing the solute, γ is the damping constant

(or collision frequency), and R(t) is the ‘white noise’ vector with mean, R(t),

of zero.

From classic theories of Brownian motion, it can be shown that although molec-

ular collisions are random, the ensemble of these collisions produces a systematic

effect. In other words, random motions exist at thermal equilibrium as a ﬂuc-

tuation. It follows that the frictional force and the random force are related by

the ﬂuctuation/dissipation theorem [684]. This relation can be expressed by the

γ-dependence of the covariance of R:

R(t)R(t



)

 =2γk

TMδ(t − t



) (10.55)

is Boltzmann’s constant and T is the temperature). Here the covariance ma-

trix is diagonal since hydrodynamic interactions between particles have been

discounted. The damping constant γ controls both the magnitude of the fric-

tional force and the variance of the random forces. It thus ensures that the system

converges to a Boltzmann distribution characterized by the temperature T. The

larger the value of γ, the greater the inﬂuence of the surrounding ﬂuctuating force

(solvent).

In the limit of small γ, the motion is termed inertial, and in the limit of

large γ it is diffusive or Brownian;see[1036] for vivid illustrations of those

regimes for models of supercoiled DNA. Different integration algorithms are gen-

erally applied depending on the relevant regime (e.g., [180,1200,1293]). See also

Chapter 14 for a simple discretization [180] of the above Langevin equation and

other discretizations in [1200,1293].

The stochastically modeled system reaches the same equilibrium as the original

system obeying M (d

x/dt

)=−∇E(x),buttherate at which equilibrium is

reached depends on the viscous coupling to the environment. Since the number of

degrees of freedom in the Langevin model is the number corresponding to the so-

lute particles, the model is computationally much cheaper than the corresponding

all-atom representation which includes explicit solvent.

Langevin Parameters from Hydrodynamic and Other Considerations

Guidelines for choosing numerical values for γ come from hydrodynamic theory.

Stokes’ law describes how the frictional resistance of a spherical particle in solu-

tion varies linearly with its radius: the effective force magnitude is 6πηa

times

the particle’s velocity, where a

is the hydrodynamic radius of each spherical par-

ticle and η is the solvent viscosity. Stokes’ law is frequently applied to particles

10.6. Continuum Solvation 337

of molecular size. Hence, γ in the Langevin equation can be set to

γ =6πηa

/m ,

where m is the particle’s mass.

From a modeling point of view, it is also possible to choose γ so that the result-

ing translational diffusion constants D

match the the experimental values in the

diffusive limit (large γ). This follows the relationship

= k



γ,

where the summation extends over all masses m

in the system. Note that this

expression for D

reduces to

= k

T/6πηa

when Stokes’ law for setting γ is used. This is the Stokes-Einstein law of diffusion

for a Brownian particle; see homework assignment 12 for an example.

A computationally efﬁcient value for γ can also be selected from practical con-

siderations, since an optimal coupling of the system to the thermal reservoir can

accelerate conﬁgurational sampling during ﬁnite-length simulations [783, 1036,

1037]. This choice of γ implies relinquishing the idea of approximating the ‘true’

dynamics in favor of efﬁcient sampling.

The Brownian Limit

In the diffusive limit of the Langevin equation, the motion is more random or

Brownian in character. Such motion characterizes in a global sense a dense system

in which the solute collides often with the surrounding ﬂuid particles, and is thus

continuously and signiﬁcantly reoriented by the solvent molecules. Speciﬁcally,

the Brownian regime assumes that the velocity relaxation time is much more rapid

than position relaxation time.

Theories for Brownian motion, dating from Einstein’s work (circa 1905), have

been based on the generalized Langevin equation and on the Fokker-Planck equa-

tion, a partial differential equation that describes the evolution of a system in a

probabilistic sense [853].

Practical BD algorithms [106, 367, 813] are similar to Monte Carlo proce-

dures (see Chapter 12) in that the current position is perturbed by a random

displacement vector. However, this random displacement is more complicated to

formulate, as it depends on the forces and the diffusion tensor. See the brief sec-

tioninChapter14 on Brownian dynamics algorithms and illustrations at the end

of Chapter 6 for long DNA and polynucleosomes.

When very large biomolecular polymers are modeled, such as long DNA of

thousands of base pairs, a more detailed account of hydrodynamic interactions

is required to accurately represent the changes in solute forces induced by the

ﬂow of the surrounding ﬂuid particles. This is done by formulating a position-

dependent hydrodynamic tensor T related to the diffusion tensor D instead of the

γ-dependent friction term in the simple Langevin equation.

338 10. Nonbonded Computations

The Brownian methods that incorporate hydrodynamic effects via a tensor that

is conﬁguration dependent yield better descriptions of large polymers in solution

and can reach longer time frames such as milliseconds [106, 107, 608]. Still, the

computational complexity grows rapidly, as O(N

), with the number of modeled

beads (N particles) when hydrodynamic effects are included [607].

Various approximations can be used to make such BD simulations with hy-

drodynamics applicable to long-time macroscopic models of biological polymers

like DNA [107]. An idea of Marshall Fixman [404] to use Chebyshev polyno-

mial approximations for matrix/vector products was developed and applied in

[1119], demonstrating that an O(N

) complexity is possible for large systems

(see Chapter 14).

10.6.4 Continuum Electrostatics

For introductions into the study of polyelectrolyte solutions and related theories,

such as the Poisson-Boltzmann equation and Debye-H¨uckel theory, readers may

wish to consult the book [537] for a good treatment of statistical thermodynamics

and kinetics, the text [547] for statistical thermodynamics, and the volume [1048]

for an early overview of how Poisson-Boltzmann and Debye-H¨uckel theories

have been applied to polyelectrolyte solutions. Classic references on electrolyte

solutions are [517, 1056]. Applications of the Poisson-Boltzmann equation for

analysis of biomolecular electrostatics are reviewed in [412] for example.

Gauss’ Law for the Electrostatic Potential

Continuum electrostatic approximations are based on numerical solutions to the

Poisson-Boltzmann (PB) equation (see [537,547], for example, for introductions).

This second-order differential equation combines theories from statistical me-

chanics — the Boltzmann distribution for a charge density ρ — with an equation

from electrostatics — Gauss’ law (or the Poisson equation),

which relates the

second derivative of the electrostatic potential Φ to the charge density.

Gauss’ law describes the electrostatic potential Φ at position x in terms of

the ﬁxed charged density of the solute, ρ

solute

(x) and the position-dependent

dielectric function (x) as:

∇·[(x) ∇Φ(x)] = −4πρ

solute

(x) . (10.56)

(See Box 10.2 for the deﬁnition of the divergence operator ∇ used above).

This equation reduces to Coulomb’s law when the dielectric constant is uniform

throughout space and the charges are modeled as point charges.

For a polar solvent like water, the effective dielectric constant increases as the

point x moves farther away from the solute. This spatially dependent dielectric

function — low  near the solute and high  for the bulk solvent — allows a better

Physicists often refer to this equation as Gauss’ law while mathematicians tend to favor the term

Poisson’s equation.

10.6. Continuum Solvation 339

description than Coulomb’s law: it takes into account the difference in electric

polarizability between the macromolecule and the solvent. (See [1420], for ex-

ample, for estimates of local dielectric constants in the environment of B-DNA in

solution). Poisson’s equation cannot be solved analytically for arbitrary geome-

tries and must be solved numerically by ﬁnite-difference or other methods (see

below).

When mobile ions are also present in the solution, the charge density is

delocalized from the solute/solvent boundary. The charge atmosphere is thus

position-dependent, since the ions redistribute in the solution in response to the

electric potential. A better approximation is obtained by considering a charge

density around the solute resulting from the distribution of charges in the medium.

Assume that we have an electrolyte solution occupying a volume V and con-

taining N

ions of corresponding charge q

for n

ion species i.Letc

= N

denote the bulk concentration of the ionic species i. The total charge density re-

sulting from the sum of all charge densities of the ions can be described by a

Boltzmann distribution as:

ρ(x)=



i=1

exp[−



(x)/k

T] , (10.57)

where



(x) is the ‘effective potential of mean force’ for electrolytes of type i at

position x for a given solute conﬁguration, or the energy change required to bring

the ion from inﬁnity to the position x.

In practice, it is assumed that



(x) is approximated by the product charge

times the potential, that is, the distribution of ions is determined by the

electrostatic ﬁeld:



(x) ≈ q

Φ(x) .

Thus, the Boltzmann distribution leads to an exponential relation between the

charge density ρ and the electrostatic potential.

The Poisson-Boltzmann Equation

Combining Gauss’ law (eq. (10.56)) with the Boltzmann charge density described

by eq. (10.57) yields the Poisson-Boltzmann (PB) equation:

∇·[(x) ∇Φ(x)] = −4πρ

solute

(x) −4π



i=1

exp[−q

Φ(x)/k

T] .

(10.58)

The PB equation is the basis for modern electrolyte theory. The solution of this

nonlinear equation for Φ(x) yields thermodynamic properties in the electrolyte

solution. Since the electrostatic potential does not vary linearly in space with

the source charges, properties of linear systems do not generally apply. As for

Poisson’s equation, analytic solutions are not available in general and various

numerical methods are used in practice.