Schlick T. Molecular Modeling and Simulation: An Interdisciplinary Guide

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

616 Appendix D. Homework Assignments

Assignment 13: Advanced Exercises in Monte Carlo

and Minimization Techniques

1. Study the function:

E(x, y)=ax

+ by

+ c (1 −cos γx)+d (1 −cos δy) . (D.11)

Note that is has many local minima and a global minimum at (x, y)=

(0, 0). Minimize E(x, y) with a =1,b=2,c=0.3, γ =3π, d =0.4,and

δ =4π by the standard simulated annealing method. Use the starting point

(1, 1) and step perturbations Δx =0.15,andsetβ in the range of 3.5 to

4.5. Limit the number of steps to ∼150. Now implement the variant of the

simulated annealing method where acceptance probabilities for steps with

ΔE<0 are proportional to exp(−βE

ΔE), with the exponent g = −1.

Analyze and compare the efﬁciency of the searches in both cases. It will be

useful to plot all pairs of points (x, y) that are generated by the method and

distinguish ‘accepted’ from ‘rejected’ points.

2. Devise a different variant of the basic simulated annealing minimization

method that would incorporate gradient information to make the searches

more efﬁcient.

3. Consider the following global optimization deterministic approach based

on the diffusion equation as ﬁrst suggested by Scheraga and colleagues

(L. Piela, J. Kostrowicki, and H. A. Scheraga, “The Multiple-Minima Prob-

lem in Conformational Analysis of Molecules. Deformation of the Potential

Energy Hypersurface by the Diffusion Equation Method”, J. Chem. Phys.

93, 3339–3346 (1989)).

The basic idea is to deform the energy surface smoothly. That is, we seek

to make “shallow” wells in the potential energy landscape disappear itera-

tively until we reach a global minimum of the deformed function. Then we

“backtrack” by successive minimization from the global minimum of the

transformed surface in the hope of reaching the global minimum of the real

potential energy surface. This idea can be implemented by using the heat

equation where T represents the temperature distribution in space x,andt

represents time:

∂

∂x

∂T

∂t

(D.12)

T (x, 0) = E(x) . (D.13)

Here, the boundary condition at time t =0equates the initial temperature

distribution with the potential energy function E(x). Under certain condi-

tions (e.g., E is bounded), a solution exists. Physically, the application of

this equation exploits the fact that the heat ﬂow (temperature distribution)

should eventually settle down.

Appendix D. Homework Assignments 617

To formulate this idea, let us for simplicity consider ﬁrst a one-dimensional

problem where the energy function E depends on a scalar x.LetE

(n)

(x)

represent the nth derivative of E with respect to x and deﬁne the

transformation operator S on the energy function E for β>0 as follows:

S[E(x)] = E(x)+βE

(2)

(x) . (D.14)

That is, we have:

E = E

E = E + βE

(2)

E = E +2βE

(2)

+ β

(4)

E = E +3βE

(2)

+3β

(4)

+ β

(8)

E =(1+βd

/dx

)

Now writing β = t/N where t is the time variable, and letting N →∞,

we write:

exp (td

/dx

) E ≡ exp(A(t)) E =



1+A +

+ ···



(D.15)

Thus we can deﬁne T (t) as

T (t)=exp(A(t)) = exp(td

/dx

) . (D.16)

In higher dimensions, let x represent the collective vector of n indepen-

dent variables; we replace the differential operator above d

/dx

by the

Laplacian operator,thatis

Δ=



i=1

∂

/∂x

Using this deﬁnition, we can also write

T (t)=T

(t) T

(t) ··· T

(t)

where

=exp(t∂

/∂x

) .

This deﬁnition produces the heat equation (D.12, D.13)since

∂T(t)[E(x)]

∂t



+ ···



[E]



1+A +

+ ···



[E]

∂

T (t)

∂x

[E(x)] .

618 Appendix D. Homework Assignments

In practice, the diffusion equation method for global optimization is

implemented by solving the heat equation by Fourier techniques (easy, for

example, if we have dihedral potentials only) or by solving for T up to

a sufﬁciently large time t. This solution, or approximate solution (repre-

senting E(x, t) for some large t), is expected to yield a deformed surface

with one (global) minimum. With a local minimization algorithm, we com-

pute the global minimum x

∗

of the deformed surface, and then begin an

iterative deformation/minimization procedure from x

∗

and E(x, t) so that

at each step we deform backwards the potential energy surface and obtain

its associated global minimum (E(x, t) −→ E(x, t − Δt) and x

∗

to x

E(x, t − Δt) −→ E(x, t − 2Δt) and x

to x

, ···E(x, 0) −→ x

). Of

course, depending on how the backtracking is performed, different ﬁnal

solutions can be obtained.

(a) To experiment with this interesting diffusion-equation approach for

global minimization, derive a general form for the deformation oper-

ator T (t)=exp(td

/dx

) on the following special functions E(x):

(i) polynomial functions of degree n,and(ii) trigonometric functions

sinωx and cosωx,whereω is a real-valued number (frequency). What

is the signiﬁcance of your result for (ii)?

(b) Apply the deformation operator T (t) = exp(td

/dx

) to the quadratic

function

E(x)=x

+ ax

+ bx

, (D.17)

with a =3and b =1. Evaluate and plot your resulting T (t)E(x)

function at t =0, Δt, 2Δt, ..., for small time increments Δt until the

global minimum is obtained.

in eq. (D.11). Examine behavior of the deformation as t →∞as a

function of the constants a and b. Under what conditions will a unique

minimum be obtained as t →∞?

4. Use Newton minimization to ﬁnd the minimum of the two-variable func-

tion in equation (D.11) and the one-variable function in equation (D.17).

It is sufﬁcient for the line search to use simple bisection: λ =1, 0.5,

etc., or some other simple backtracking strategy. For the quartic function,

experiment with various starting points.

Read remaining paper from Coursepack

Appendix D. Homework Assignments 619

Assignment 14: Advanced Exercises in Molecular

Dynamics

1. Calculate the ‘natural’ time unit for molecular dynamics simulations of

biomolecules from the relation: energy = mass ∗ (length/time)

, to obtain

the time unit τ corresponding to the following units:

length (l): 1

A =10

−10

mass (m): 1 amu = 1g/mol

energy (v): 1 kcal/mol = 4.184 kJ/mol.

Estimate the “quantum mechanical cutoff frequency”, ω

= kT/ at room

temperature (∼300

◦

K).

2. Derive the amplitude decay rate of γ/2 for an underdamped harmonic

oscillator due to friction by solving the equations of motion:

= −kx − mγ

(D.18)

and examining time behavior of the solution.

3. Derive the amplitude decay rate of ω

(Δt)/2 intrinsic to the implicit-Euler

scheme by solving the discretized form of eq. (D.18).

4. Compare your answer in problem 2 above with behavior of the explicit-

Euler solution of eq. (D.18).

5. Compare molecular and Langevin dynamics simulations of two water mole-

cules by the Verlet discretization of the equation of motion and its Langevin

analog. Use the “SPC” intermolecular potential, given by:

E =



(i,j)≡(O,O) pairs

i<j



−A





(k,)≡intermolecular

(O,O),(O,H),(H,H) pairs

k<





k



where

A = 626 (kcal

)/mol

B = 629 × 10

(kcal

)/mol

=0.41 e

= −0.82 e.

A numerical factor of 332 is needed in the electrostatic potential to ob-

tain energies in kcal/mol with the coefﬁcients above. For simplicity, assume

that intramolecular geometries are rigid: r

= 1

A, cos θ

HOH

= −1/3.

(You can use harmonic soft constraints). Begin by ﬁrst minimizing the

energy of the water dimer and examining the hydrogen bond geometry

620 Appendix D. Homework Assignments

(hydrogen-bond distance and angle θ between the hydrogen-bond vector

and bisector of the acceptor molecule). Then study numerical behavior of

the two models/schemes as a function of Δt, and examine the hydrogen

bond geometry. Experiment with Δt =1,2,5,and10fsanduseγ val-

ues in the range of 1 to 50 ps

−1

. If you are more ambitious, continue to

study energy-minimized structures of water clusters of larger sizes and their

dynamics. Analyze the hydrogen bonding networks of these clusters.

Donor Acceptor

Figure D.4. Hydrogen bond geometry: the angle θ is deﬁned between the hydrogen-bond

vector and the bisector of the acceptor molecule.

Some Useful Constants and Conversion Factors

Avogadro’s Number N

=6.0221 × 10

mol

−1

Planck’s Constant h =6.6261 × 10

−34

Jsec

 = h/2π =1.055 × 10

−34

Jsec

Boltzmann’s Constant k

=1.38066 × 10

−23

−1

Gas Constant R = k

=8.3145 JK

−1

mol

−1

Atomic Mass Unit, amu (1/N

)=1g/mol = 1.6605 ×10

−27

π =3.14159

1 kcal = 4.184 kJ

Appendix D. Homework Assignments 621

Assignment 15: BONUS PROJECT

The Scaling of Protein-Conformer Number with Size and Solvability of the

Protein Folding Problem

The phone rings one morning as you sip through your Sunluck vanilla latt´e

grande and check your emails in your ofﬁce at the University of Seabeetle. The

editor of the journal Proteomics Today is on the line.

Given your expertise in biomolecular modeling, she asks your help in writing

abriefFolding In Silico Perspectives article for the next issue discussing the fol-

lowing series of interesting papers debating the nature of scaling of the number

of protein conformers as function of chain length (exponential or nonexponential)

and the solvability of the protein folding problem by computer simulation:

• W. F. van Gunsteren, R. B¨urgi, C. Peter, and X. Daura, “The Key to Solv-

ing the Protein-Folding Problem Lies in an Accurate Description of the

Denatured State”, Angew. Chem. Int. Ed. 40, 352–355 (2001).

• A. R. Dinner and M. Karplus, “Comment on the Communication ‘The Key

to Solving the Protein-Folding Problem Lies in an Accurate Description of

the Denatured State’ by van Gunsteren et al.”, Angew. Chem. Int. Ed. 40,

4615–4616 (2001).

• W. F. van Gunsteren, R. B¨urgi, C. Peter, and X. Daura, “Reply”, Angew.

Chem.Int.Ed.40, 4616–4618 (2001).

Since many of Proteomics Today readers and authors perform computer sim-

ulations of proteins, both macroscopic and all-atom based, the editor wants you

also to mention the strengths and weaknesses of these different approaches.

Though already overloaded with preparing ﬁnal examinations for your classes,

writing grant proposals, supervising your students and postdocs, and reviewing

several articles for journals (long-overdue), you agree to take this challenging and

potentially rewarding assignment. You immerse yourself in the papers, send your

graduate students to collect background articles and information, order through

your cellphone Sunluck’s caramel-macchiato enormoso, and decide to make your

article not only objective and interesting but also fun to read.

Some Suggestions:

• Discuss the Levinthal paradox.

• What is the relation among the number of conformers, timescales, and

folding pathways?

• Analyze lattice or other protein simulations reported in the literature to

estimate the number of possible conformers and the timescale of protein

folding.

References

[1] L. Adams and J. L. Nazareth, editors. Linear and Nonlinear Conjugate Gradient-

Related Methods. SIAM, Philadelphia, PA, 1996.

[2] M. D. Adams, G. G. Sutton, H. O. Smith, E. W. Myers, and J. C. Venter.

The independence of our genome assembly. Proc. Natl. Acad. Sci. USA, 100:

3025–3026, 2003.

[3] P. L. Adams, M. R. Stahley, A. B. Kosek, J. Wang, and S. A. Strobel. Crystal

structure of a self-splicing group I intron with both exons. Nature, 430:45–50,

2005.

[4] S. A. Adcock and J. A. McCammon. Molecular dynamics: survey of methods for

simulating the activity of proteins. Chem. Rev., 106:1589–1615, 2006.

[5] D. K. Agraﬁotis. Stochastic algorithms for maximizing molecular diversity.

J. Chem. Inf. Comput. Sci., 37:841–851, 1997.

[6] D. K. Agraﬁotis. Diversity of chemical libraries. In P. von Ragu´e Schleyer (Editor-

in Chief), N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, and H. F. Schaefer,

III, editors, Encyclopedia of Computational Chemistry, volume 1, pages 742–761.

John Wiley & Sons, West Sussex, England, 1998.

[7] D. K. Agraﬁotis, V. S. Lobanov, and F. R. Salemme. Combinatorial informatics in

the post-genomics era. Nat. Rev. Drug. Disc., 1:337–346, 2002.

[8] D. K. Agraﬁotis, J. C. Myslik, and F. R. Salemme. Advances in diversity proﬁling

and combinatorial series design. Mol. Div., 4:1–22, 1999.

[9] A. Aguzzi. Prions and antiprions. Biol. Chem., 378:1393–1395, 1997.

[10] A. Aguzzi, F. Montrasio, and P. S. Kaeser. Prions: Health scare and biological

challenge. Nature Rev. Mol. Cell Biol., 2:118–126, 2001.

[11] P. Ahlquist. RNA-dependent RNA polymerases, viruses, and RNA silencing.

Science, 296:1270–1273, 2002.

624 References

[12] J. Aishima, R. K. Gitti, J. E. Noah, H. H. Gan, T. Schlick, and C. Wolberger.

A Hoogsteen base pair embedded in undistorted B-DNA. Nuc. Acids Res.,

30:5244–5252, 2002.

[13] J. Aishima and C. Wolberger. Crystal structure of the MATα2 homeodomain-DNA

complex with nonspeciﬁcally bound homeodomains. Nuc. Acids Res., 2002.

[14] A. Aiuti, F. Cattaneo, S. Galimberti, U. Benninghoff, B. Cassani, L. Callegaro,

S. Scaramuzza, G. Andolﬁ, M. Mirolo, I. Brigida, A. Tabucchi, F. Carlucci,

M. Eibl, M. Aker, S. Slavin, H. Al-Mousa, A. Al Ghonaium, A. Ferster,

A. Duppenthaler, L. Notarangelo, U. Wintergerst, R. Buckley, M. Bregni,

S. Marktel, M. Valsecchi, P. Rossi, F. Ciceri, R. Miniero, C. Bordignon, and

M. Roncarolo. Gene therapy for immunodeﬁciency due to adenosine deaminase

deﬁciency. New Engl. J. Med., 360:447–458, 2009.

[15] E. Akhmatskaya, N. Bou-Rabee, and S. Reich. A comparison of generalized hy-

brid Monte Carlo methods with and without momentum ﬂip. J. Comput. Phys.,

228:2256–2265, 2009.

[16] B. Al-Lazikani, J. Jung, Z. Xiang, and B. Honig. Protein structure prediction. Curr.

Opin. Struct. Biol., 5:51–56, 2001.

[17] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter. Molecu-

lar Biology of the Cell. Garland Science Publishing, New York, NY, ﬁfth edition,

2008.

[18] I. L. Alberts, Y. Wang, and T. Schlick. DNA polymerase β catalysis: Are different

mechanisms possible? J. Amer. Chem. Soc., 129:11100–11110, 2007.

[19] B. J. Alder and T. E. Wainwright. Studies in molecular dynamics. I. General

method. J. Chem. Phys., 31:459–466, 1959.

[20] C. Alhambra, F. J. Luque, F. Gago, and M. Orozco. Ab Initio study of stacking

interactions in A- and B-DNA. J. Phys. Chem. B, 101:3846–3853, 1997.

[21] J. F. Allemand, D. Bensimon, R. Lavery, and V. Croquette. Stretched and over-

wound DNA forms a Pauling-like structure with exposed bases. Proc. Natl. Acad.

Sci. USA, 95:14152–14157, 1998.

[22] M. P. Allen and D. J. Tildesley. Computer Simulation of Liquids. Oxford University

Press, New York, NY, 1990.

[23] R. J. Allen, D. Frenkel, and P. R. ten Wolde. Simulating rare events in equilibrium

or nonequilibrium stochastic systems. J. Chem. Phys., 124:024102, 2006.

[24] N. L. Allinger. Calculation of molecular structure and energy by force-ﬁeld

methods. Adv. Phys. Org. Chem., 13:1–85, 1976.

[25] N. L. Allinger. Conformational analysis. 130. MM2 A hydrocarbon force ﬁeld

utilizing V1 and V2 torsional terms. J. Amer. Chem. Soc., 99:8127–8134, 1977.

[26] N. L. Allinger and K. Chen. Hyperconjugative effects on carbon-carbon bond

lengths in molecular mechanics (MM4). J. Comput. Chem., 17:747–755, 1996.

[27] N. L. Allinger, J. T. Fermann, W. D. Allen, and H. F. Schaefer, III. The torsional

conformations of butane: Deﬁnitive energetics from ab initio methods. J. Chem.

Phys., 106:5143–5150, 1997.

[28] N. L. Allinger, M. A. Miller, F. A. VanCatledge, and J. A. Hirsch. Conformational

analysis. LVII. The calculation of the conformational structures of hydrocar-

bons by the Westheimer-Hendrickson-Wiberg method. J. Amer. Chem. Soc.

, 89:

4345–4357, 1967.

References 625

[29] N. L. Allinger and J. T. Sprague. Calculation of the structures of hydrocarbons

containing delocalized electronic systems by the molecular mechanics method.

J. Amer. Chem. Soc., 95:3893–3907, 1973.

[30] N. L. Allinger, Y. H. Yuh, and J.-H. Lii. Molecular mechanics. The MM3 force ﬁeld

for hydrocarbons. 1. J. Amer. Chem. Soc., 111:8551–8566, 1989.

[31] C. D. Allis, T. Jenuwein, and D. Reinberg, editors. Epigenetics. Cold Spring

Harbor Laboratory Press, Cold Spring Harbor, NY, 2009.

[32] S. A. Allison and J. A. McCammon. Multistep Brownian dynamics: Application

to short wormlike chains. Biopolymers, 23:363–375, 1984.

[33] U. Alon. Network motifs: theory and experimental approaches. Nat. Rev. Genet.,

8:450–461, 2007.

[34] A. Altis, M. Otten, P. H. Nguyen, R. Hegger, and G. Stock. Construction of the free

energy landscape of biomolecules via dihedral angle principal component analysis.

J. Chem. Phys., 128:245102, 2008.

[35] C. Altona and M. Sundaralingam. Conformational analysis of the sugar ring in nu-

cleosides and nucleotides. a new description using the concept of pseudorotation.

J. Amer. Chem. Soc., 94:8205–8212, 1972.

[36] A. Amadei, A. B. M. Linssen, and H. J. C. Berendsen. Essential dynamics of

proteins. Proteins: Struc. Func. Gen., 17:412–425, 1993.

[37] A.Amadei,A.B.M.Linssen,B.L.deGroot,D.M.F.vanAalten,and

H. J. C. Berendsen. An efﬁcient method for sampling the essential subspace of

proteins. J. Biomol. Struct. Dynam., 13:615–625, 1996.

[38] J. G. Amar. The Monte Carlo method in science and engineering. Comput. Sci.

Engin., 8:9–19, 2006. doi:10.1109/MCSE.2006.34.

[39] A. Amir-Aslani. Toxicogenomic predictive modeling: Emerging opportunities for

more efﬁcient drug discovery and development. Tech. Forecast. Soc. Change,

75:905–932, 2008.

[40] M. Amos. Theoretical and Experimental DNA Computation. Natural Computing

Series. Springer, New York, NY, 2005.

[41] L. M. Amzel. Structure-based drug design. Curr. Opin. Biotech., 9:366–369, 1998.

[42] K. Anand, J. Ziebuhr, P. Wadhwani, J. R. Mesters, and R. Hilgenfeld. Coronavirus

main proteinase (3CL

pro

) structure: Basis for design of anti-SARS drugs. Science,

300:1763–1767, 2003.

[43] H. C. Andersen. Molecular dynamics simulations at constant pressure and/or

temperature. J. Chem. Phys., 72:2384–2393, 1980.

[44] H. C. Andersen. Rattle: a ‘velocity’ version of the SHAKE algorithm for molecular

dynamics calculations. J. Comput. Phys., 52:24–34, 1983.

[45] C. R. Anderson. An implementation of the fast multipole algorithm without

multipoles. SIAM J. Sci. Stat. Comput., 13:923–947, 1992.

[46] S. L. Anderson. Random number generators on vector supercomputers and other

advanced architectures. SIAM Rev., 32:221–251, 1990.

[47] A. Andreeva, D. Howorth, S. E. Brenner, T. J. P. Hubbard, C. Chothia, and

A. G. Murzin. SCOP database in 2004: Reﬁnements integrate structure and

sequence family data. Nucl. Acids Res., 32:D226–D229, 2004.

626 References

[48] A. Andreeva, D. Howorth, J.-M. Chandonia, S. E. Brenner, T. J. P. Hubbard,

C. Chothia, and A. G. Murzin. Data growth and its impact on the SCOP database:

New developments. Nucl. Acids Res., 36:D419–D425, 2008.

[49] G. Andres Cisneros, L. Perera, M. Garcia-Diaz, K. Bebenek, T. A Kunkel, and

L. G. Pedersen. Catalytic mechanism of human DNA polymerase λ with Mg

and Mn

from ab initio quantum mechanical/molecular mechanical studies. DNA

Repair, 7:1824–1834, 2008.

[50] M. Andronescu, V. Bereg, H.H. Hoos, and A. Condon. RNA STRAND: the RNA

secondary structure and statistical analysis database. BMC Bioinformatics, 9:340,

2008.

[51] C. B. Anﬁnsen, E. Haber, M. Sela, and F. H. White, Jr. The kinetics of formation of

native ribonuclease during oxidation of the reduced polypeptide chain. Proc. Natl.

Acad. Sci. USA, 47:1309–1314, 1961.

[52] C. B. Anﬁnsen and H. A. Scheraga. Experimental and theoretical aspects of protein

folding. Adv. Protein Chem., 29:205–301, 1975.

[53] A. W. Appel. An efﬁcient program for many-body simulation. SIAM J. Sci. Stat.

Comput., 6:85–103, 1985.

[54] J.

Aqvist and A. Warshel. Simulation of enzyme reactions using valence bond force

ﬁelds and other hybrid quantum/classical approaches. Chem. Rev., 93:2523–2544,

1993.

[55] G. Arents, R. W. Burlingame, B. C. Wang, W. E. Love, and E. N. Moudrianakis.

The nucleosomal core histone octamer at 3.1

A resolution: A tripartite protein as-

sembly and a left-handed superhelix. Proc. Natl. Acad. Sci. USA, 88:10148–10152,

1991.

[56] A. Arkhipov, P. L. Freddolino, and K. Schulten. Stability and dynamics of virus

capsids described by coarse-grained modeling. Structure, 14:1767–1777, 2006.

[57] A. Arkhipov, Y. Yin, and K. Schulten. Four-scale description of membrane

sculpting by BAR domains. Biophys. J., 95:2806–2821, 2008.

[58] Y. A. Arnautova, A. Jagielska, and H. A. Scheraga. A new force ﬁeld (ECEPP-05)

for peptides, proteins, and organic molecules. J. Phys. Chem. B, 110:5025 W5044,

2006.

[59] K. Arora and T. Schlick. Deoxyadenosine sugar puckering pathway simulated by

the stochastic difference equation algorithm. Chem. Phys. Lett., 378:1–8, 2003.

[60] K. Arora and T. Schlick. In Silico evidence for DNA polymerase β’s substrate-

induced conformational change. Biophys. J., 87:3088–3099, 2004.

[61] P. G. Arscott, G. Lee, V. A. Bloomﬁeld, and D. F. Evans. Scanning tunneling

microscopy of Z-DNA. Nature, 339:484–486, 1989.

[62] G. Arya and T. Schlick. Role of histone tails in chromatin folding revealed

by a mesoscopic oligonucleosome model. Proc. Natl. Acad. Sci. USA,

103:16236–16241, 2006.

[63] G. Arya and T. Schlick. Efﬁcient global biopolymer sampling with end-transfer

conﬁgurational bias Monte Carlo. J. Chem. Phys., 106:044107, 2007.

[64] G. Arya and T. Schlick. A tale of tails: How histone tails mediate chromatin

compaction in different salt and linker histone environments. J. Phys. Chem. A,

113:4045–4059, 2009.