Molekulová mechanika empirické potenciály silová pole
Michal Otyepka, PřF UP Olomouc
Proč, když máme QM? běžná malá molekula – kvantový chemik jásá středně velká molekula – kvantovému chemikovi tuhnou rysy a volí umírněné prostředky (často i HF/6-31G(d))
běžný protein – kvantový chemik je v pasti
Born-Oppenheimerova aproximace oddělení elektronického a jaderného pohybu kvantové elektrony vs. klasická jádra
E = f (R ) aparát klasické fyziky molekulová mechanika
1
Potenciální energie
Deformace vazebné vzdálenosti zajímá nás chování okolo minima využijeme Taylorův rozvoj -10
vzdálenost 10 1.0
m
2.0
elektronická energie
( )
E (r ) = E ro +
0.74 = r0 vazebná vzdálenost
(
( ) (r − r ) + 1 ∂ E (r ) (r − r ) 2!
1 ∂E r0 1! ∂r
k E ( r ) = r − ro 2
2
2
0
0
0
∂r 2
+ ...
)
2
Vazba jako pružina ∂E (r ) síla = −k r − r0 ∂r výchylka z rovnovážné polohy r − r0 = x
(
E
F =−
F = −kx E (r ) =
(
k r − ro 2
)
Hookův zákon silová konstanta
)
2
∂ 2 E (r ) ∂r 2
=k
tato aproximace nedovolí disociaci vazeb!
r
r0
E=
(
)
2 1 ED 1 − e −a ( r − r0 ) 2 Morseho pot.
Různé vazby různé pružiny různé kovalentní vazby mají různou vazebnou vzdálenost silovou konstantu molekula H2 H35Cl H79Br H127I
k / N·m-1 510
r0 / pm 74.1
478 408 291
127.5 141.4 160.9
2
Typy vazeb podobné vazby X–Y se chovají ve všech molekulách podobně bez ohledu na okolí - parametry jsou přenositelné najít podobné vazby, přiřadit jim k a r0 zavádí se atomové typy
Atomové typy - uhlík PARM99 for DNA,RNA,AA, organic molecules, TIP3P wat. C sp2 C carbonyl group CA sp2 C pure aromatic (benzene) CB sp2 aromatic C, 5&6 membered ring junction CC sp2 aromatic C, 5 memb. ring HIS CD sp2 C atom in the middle of: C=CD-CD=C CK sp2 C 5 memb.ring in purines CM sp2 C pyrimidines in pos. 5 & 6 O CN sp2 C aromatic 5&6 memb.ring junct.(TRP) CQ sp2 C in 5 mem.ring of purines between 2 N CR sp2 arom as CQ but in HIS CT sp3 aliphatic C CV sp2 arom. 5 memb.ring w/1 N and 1 H (HIS) CW sp2 arom. 5 memb.ring w/1 N-H and 1 H (HIS) C* sp2 arom. 5 memb.ring w/1 subst. (TRP) CY nitrile C (Howard et al.JCC,16,243,1995) CZ sp C (Howard et al.JCC,16,243,1995)
H HC C
CT CT
O
H
HC
H HC N H
N
H
Ala
Vazebné typy
CT-CT CT-HC CT-H1 CT-H2 CT-H3 CT-HP CT-N* CT-N2 CT-OH CT-OS C*-HC C*-CB C*-CT C*-CW CB-CN
k
r0
310.0 340.0 340.0 340.0 340.0 340.0 337.0 337.0 320.0 320.0 367.0 388.0 317.0 546.0 447.0
1.526 1.090 1.090 1.090 1.090 1.090 1.475 1.463 1.410 1.410 1.080 1.459 1.495 1.352 1.419
databáze parametrů – pole pole balíku AMBER
3
Jak získat parametry? z experimentů
vazebné geometrie z
vibrační spektroskopie
výpočtem
IR data (free) NIST - http://webbook.nist.gov/chemistry/
RTG a neutronová difrakce, NMR, rotační spektroskopie
silové konstanty z
Struktury (free) Protein Data Bank (PDB) – http://www.pdb.org
ν~ =
1 ⎛⎜ k 2πc ⎜⎝ meff
1/ 2
⎞ ⎟ ⎟ ⎠
meff =
m1 m2 m1 + m2
fitováním energetických hyperploch vypočtených referenční QM metodou
Molekulová mechanika celková energie je funkcí vzájemné pozice jader
E = f (R ) = Ecovalent + Enoncovalent
Ecovalent = Eb + Ea + Et Enoncovalent = Ec + Evdw aditivní model
Deformace úhlů E=
(
)
2 kθ θ − θ0 2 kθ = 80 kcal/mol. deg 2
θ 0 = 122.9°
4
Deformace torzí
2.9 kcal/mol
E
HH
60
0
H H HH
H
120
H
H
H H
180
240
Degrees of Rotation θ
H
300
360
E=
(
(
kt 1 + cos nφ − φ0 2
))
n=3 φ0 = 180.0° kt = 2.9/2*9 (IDIVF=1) kt = 2.9/2 (IDIVF=9)
E=∑ i
(
(
k i ,t 1 + cos ni φi − φi , 0 2
))
Konformační chování konformační chování biomakromolekul je do značné míry dáno pamaterizací torzní úhlů (a nevazebnými termy) parametrizace
podle experimentálních dat (NMR, distribuce z databází) – „knowledge based“ – CHARMM QM profily – QM v gas phase není zcela vhodná, QM v solventu (εr = 80) je lepší alternativa
CMAP – 2D dihedral energy correction map, CHARMM
Deformace torzí
C-C-C-C H-C-C-H
5
Nepřímé torze např. pro popis vhodné geometrie aminoskupin
AMBER 2-3-1-4 fáze 180° n=2 CHARMM
E=
(
kω ω − ω0 2
)
2
Nekovalentní interakce (
)
(
) ∑ u (R , R , R )+ ... ) model párového potenciálu
U R1 ,..., R N = ∑ u R i , R j +
(
)
i< j
(
U R1 ,..., R N = ∑ u R i , R j
(
)
i< j
(
)
i
i< j
(
u R i , R j = u atr R i , R j + u rep R i , R j
j
k
)
Elektrostatická interakce popisuje interakci multipól-multipól v monopólovém rozvoji atomově centrované parciální náboje Coulombův zákon
E=
1
qi q j
4πε 0
r
6
Elektrostatika
Parciální náboje Mullikenovy – nevhodné RESP náboje – „Restrained ElectroStatic Potential fit“ – dobře fitují elektrostatický potenciál molekuly
HF/6-31G* - přecenění dip. momentů, částečná kompenzace elektrostatické indukce (parm99) B3LYP/cc-pVTZ/PCM(εr=4) … (ff03) polarizaci lze zavést dodatečně
RESP
7
RESP
Van der Waalsova interakce popisuje disperzní a repulzní složku nekovalentní interakce výhodný je popis Lennard-Jonesovým potenciálem
Energie (kJ/mol)
0.2
0.1
repulze
0
hloubka minima
vdW profil
-0.1
disperze -0.2 2
3
4
5
6
7
8
9
10
vzdálenost (10-10 m)
Lennard-Jonesův potenciál ⎛ ⎛ σ ⎞12 ⎛ σ ⎞ 6 ⎞ ⎛ ⎛ σ ⎞12 ⎛ σ ⎞ 6 ⎞ u (r ) = 4ε ⎜ ⎜ ⎟ − ⎜ ⎟ ⎟ = ε ⎜ ⎜⎜ vdw ⎟⎟ − 2⎜⎜ vdw ⎟⎟ ⎟ ⎜⎝ r ⎠ ⎜⎝ r ⎠ ⎝ r ⎠ ⎟⎠ ⎝ r ⎠ ⎟⎠ ⎝ ⎝
(
)
u (σ ) = 0, u 6 2σ = min = ε ,
σ vdw / 2 van der Waalsův poloměr
0.2
Proč je výhodný LJ potenciál 12-6? Počítač umí rychle počítat mocniny a r –12 = (r –6)2.
f(r –12) 0.1
0
-0.1
f(r –6) -0.2 2
3
4
5
6
7
8
9
10
8
Molekulová mechanika (
)
(
)
2 kr r − r0 2 2 k Ea = θ θ − θ 0 2
Eb =
(
(
kt 1 + cos nφ − φ0 2 1 qi q j Ec = 4πε 0 ε r rij
Et =
⎛σ Evdw = −2ε ij ⎜ ij ⎜r ⎝ ij
))
6
⎛σ ⎞ ⎟ + ε ij ⎜ ij ⎜r ⎟ ⎝ ij ⎠
12
⎞ ⎟ ⎟ ⎠
Kompenzace chyb – parm99 Minimum – vdw komplex
Závislost rozdílu jednotlivých příspěvků interakční energie, počítaných na úrovni MM s příspěvky, počítaných na úrovni SAPT, na vzdálenosti. Eelst elektrostatická energie, Edisp disperzní energie, Eexch repulzní energie, Etot celková energie. Pro komplex ma+...mOH a ff99.
Topologie molekuly definuje, které vazby, úhly, torze etc. se uplatňují v molekule ALA
INT
1
CORR OMIT DU
BEG
0.00000 1
DUMM
DU
M
0
-1
-2
0.000
0.000
0.000
0.00000
2
DUMM
DU
M
1
0
-1
1.449
0.000
3
DUMM
DU
M
2
1
0
1.522
111.100
4
N
N
M
3
2
1
1.335
116.600
0.000
0.00000
0.000
0.00000
5
H
H
E
4
3
2
1.010
119.800
0.000
0.27190
6
CA
CT
M
4
3
2
1.449
121.900
180.000
0.03370
7
HA
H1
E
6
4
3
1.090
109.500
300.000
0.08230
8
CB
CT
3
6
4
3
1.525
111.100
60.000 -0.18250
9
HB1
HC
E
8
6
4
1.090
109.500
60.000
0.06030
HB2
HC
E
8
6
4
1.090
109.500
180.000
0.06030
11
HB3
HC
E
8
6
4
1.090
109.500
300.000
0.06030
12
C
C
M
6
4
3
1.522
111.100
180.000
0.59730
13
O
O
E
12
6
4
1.229
120.500
název at. typ konekt.
vzdálenost úhel
12 13
180.000 -0.41570
10
H9
H7
8 6
H10
O H11 N H
5
4
Ala
0.000 -0.56790
torze
parc. náboj
9
Ne vše se počítá vazebné – jen kovalentně vázaní sousedé vaz. úhel – jen reálné vazebné úhly torze – jen reálné torze coulomb – 1-2, 1-3 se nepočítají; 1-4 se škálují (2.0), další všechny 2 vdW – 1-2 a 1-3 se nepočítají; 1-4 se 1 škálují (1.2), další všechny
3 4
snížení počtu nekov. interakcí – zavádí se cutoff
Ořezání (cutoff) počet nevazebných interakcí ~ 2 × N(N–1)/2
pro 10.000 atomů (menší systém) ~108 párů
vdW interakce velmi rychle vyhasíná se vzdáleností (r–6)
vdW interakce v 2.4násobné vzdálenosti než odpovídá minimu je cca 100× menší než v minimu zavádí se cutoff pro vdw interakci; páry se vzdáleností nad rmax se nepočítají problém: vnesení diskontinuity (může způsobovat problémy tam, kde se pracuje s derivacemi E) z
switch
Cutoff
T. Sprik
10
Cutoff pro elektrostatiku? elektrostatická interakce dvou E = tot monopólů vyhasíná pomalu (r–1)
na vzdálenosti 10r činí 10% vzhledem k interakci na vzdálenosti r
∑ φ (r
i
i, j
)
− r j = Esr + Elr
(
Esr = ∑ φsr ri − r j
)
i, j
~ (k ) ρ~ (k ) 2 při použití periodických Etr = Φ lr okrajových podmínek (PBC) k lze řešit Ewaldovou sumací ~ PME – particle mesh Ewald Φ lr Fourierův obraz potenciálu
Interakce dvou částic se sumuje v „krátkých“ vzdálenostech přímo a v „dlouhých“ vzdálenostech v Fourierově prostoru a využívá Fast Fourier Transform
∑
ρ~ (k ) Hustota náboje
Další zjednodušení popis celých skupin (united atom)
např. skupina –CH3 se bude popisovat jako jeden „pseudoatom“
Coarse grained modely (CG)
Studium velkých systémů Hrubá aproximace Aplikace: membrány, proteinové komplexy CGFF - MARTINI
Coarse Grained MARTINI
The forcefield has been parametrized in a systematic way, based on the reproduction of partitioning free energies between polar and apolar phases of a large number of chemical compounds. The model is based on a four-to-one mapping, i.e. on average four heavy atoms are represented by a single interaction center. In order to keep the model simple, only four main types of interaction sites are defined: polar (P), non-polar (N), apolar (C), and charged (Q). Each particle type has a number of subtypes, which allow for an accurate representation of the chemical nature of the underlying atomistic structure. Currently topologies are available for many lipids and surfactant molecules, including cholesterol, and for all amino acids. Scripts are furthermore available to build topologies for arbitrary peptides and proteins.
11
CG - MARTINI
CG - MARTINI
Běžná silová pole organické molekuly
MM2, MM3, CVFF ...
biomakromolekuly
AMBER z
parm94, 98, ff03
CHARMM OPLS+
12
AMBER
Molekulová dynamika klasická molekulová dynamika
Molekulová dynamika časový vývoj systému
13
Časové škály proteinových pohybů průchod iontu kanálem Porinu elastické pohyby proteinů
skládání α-šroubovice skládáníβ-hairpinu
vibrace vazeb
skládání proteinu
10-15
10-12
10-9
10-6
10-3
100
(fs)
(ps)
(ns)
(μs)
(ms)
(s)
Čas
MD
Molekulová dynamika Fi = mi a i = mi
∂ 2ri ∂t
2
Fi = −
= mi &r&i
integrace Newtonových rovnic ri (t + Δt ) = ri (t ) + (Δt )v i (t ) + nová
stará
stará
∂E ri ∂ri ri
potenciál MM
stará 1 (Δt )2 Fi (t ) 2 mi
velocity Verlet
aktualizované souřadnice v čase Δt (1-2 fs)
Aktualizace rychlostí F (t ) 1 ⎛ Δt ⎞ v i ⎜ t + ⎟ = v i (t ) + (Δt ) i 2 mi 2⎠ ⎝
rychlosti se aktualizují ve dvou krocích Δt/2
F (t + Δt ) ⎛ Δt ⎞ 1 v i (t + Δt ) = v i ⎜ t + ⎟ + (Δt ) i 2⎠ 2 mi ⎝
teplota a primitivní realizace termostatu N 1 N 1 3 mi v i2 ≅ Nk BT T= ∑ ∑ mi v i2 2 3 Nk B i =1 i =1 2
14
Trajektorie
stav 0 ri0
stav 1 ri1
stav 2
stav 3
F i = mi a i
snímky = MD trajektorie (ri,j, vi,j) díky Evo
termostat (T)
Jak začít? generace náhodných rychlostí pro danou teplotu
⎛ v2 ⎞ exp⎜ − i 2 ⎟, N 0, σ 2 ⎜ 2σ ⎟ ⎝ ⎠ k T 2 σ 2 = v i2 − v i = v i2 = B mi p( v i ) =
1
σ 2π
(
)
v systému musí být malé gradienty, jinak může dojít k „explozi“ systému – před MD simulací je třeba systém minimalizovat výhodné je také postupné zvyšování teploty
Směrem k realitě neutralita
vyžadována elektroneutralita systému (náboj solutu se kompenzuje přidáním protiiontů)
solvent
solut je obklopen solventem (nejčastěji vodou) explicitní vs. implicitní model
periodické okrajové podmínky
15
Periodické okrajové podmínky ∞
∞
v
∞
∞
Periodické okrajové podmínky ∞
∞
∞
∞
Mašinérie výpočtu struktura solutu – PDB databáze, etc. příprava solutu – protonace, dostavba chybějících částí, parametrizace nestandardních reziduí etc. solvatace a přidání protiiontů
pozn. ke krystalovým vodám proteinů
simulační protokol výpočet (AMBER, GROMACS, CHARMM ...) analýza
vizualizace VMD, gOpenMol, MoilView ...
16
Solut - protein
Solut + protiionty
Solut + protiionty + voda (box)
93.678 atomů
17
Periodické okrajové podmínky
Water molecule V-shaped tetrahedrally arranged sp3-hybridized electron pairs
Implicit water models Solvent is treated as a polarizable continuum with a dielectric constant [ε] instead of explicit solvent molecules. Significantly cheaper than explicit solvent models Cannot model specific interactions such as hydrogen bonds
18
Mixed water models The first solvation sphere is explicitly described by a number of solvent molecules. The remaining solvent molecules are described by an uniform continuum medium with a dielectric constant.
ε
Explicit models rigid models (fixed atom positions) flexible models (include bond streching & angle bending) polarizable models (include explicit polarization term) planar steric
3-SITE planar & rigid three interaction sites each atom gets assigned a point charge TIP3P (Jorgensen,1983) SPC/E (Berendsen,1987) POL3 (Caldwell & Kollman, 1995)
19
4-SITE planar & rigid the negativ charge near the oxygen along the bisector of the HOH angle TIP4P (Jorgensen,1983)
5-SITE the negative charge representing the lone pairs of the oxygen atom tetrahedral & rigid TIP5P (Jorgensen & Mahoney, 2000)
Analýza – co je RDF? Bunches of red oxygens and white hydrogens show the probability that a water molecule at certain distance r from the fixed (blue and yellow) molecule will have given relative position and orientation. The animation is thus an attempt to visualize the full pair correlation function. This function is six-dimensional: one distance r and five angular variables. We have measured this function during a 300 ps simulation of 500 TIP4P water molecules at T=25oC and ambient density. The histograms were taken by 15o in each angular variable (with the exception of regions close to the poles). The r-coordinate is identified with time, i.e., each frame corresponds to a selected distance (shown in red in the graph of the O-O pair correlation function). For this distance, one molecule from the pair is fixed, and around it 8000 positions of water molecules with highest probability are drawn with diameters proportional to the probability. See also the water page.
Kolafa
20
Radial distribution function
Radiální Distribuční Funkce
v systému solut-solvent – první solvatační sféra
strukturní faktor je spojen s Fourierovou transformací RDF
Radial distribution function – 1atm 1,9
3,5
1,7
TIP3P
1,5
3
TIP4P
TIP5P
SPC/E
POL3
1,3 1,1
2,5
0,9 0,7
2
0,5 3
3,5
4
4,5
5
5,5
6
1,5 1 0,5 0 0
2
4
6
8
10
21
Explicit water models Dipole moment
Dielectric constant
Self diffusion [10-5cm2/s]
Density maximum [ °C]
Melting temperatures [ °C]
SPC/E
2.4
71
2.5
-38.0
-58.2
TIP3P
2.4
82
5.2
-91.0
-127.6
TIP4P
2.2
53
3.3
-25.0
-40.7
TIP5P
2.3
81.5
2.6
+4.0
0.8
Expt.
2.95
78.4
2.3
+3.984
0
TIP3P byl parametrizován pro cutoff, při použití PME se zvyšuje D, reparametrizaci navrhl Brooks (J. Chem. Phys. 121 (20), 10096, 2004
Cost of computation
TIP3P SPC/E TIP4P TIP5P POL3
100% ~101% ~140% ~240% ~1230%
Simulační protokol minimalizace přidaných vodíků solutu minimalizace protiiontů a vod krátká NpT simulace protiiontů a vod, dokud g=1 g/cm3 minimalizace solutu termalizace, NpT simulace s pomalu rostoucí teplotou např. k 298.15K (300K) vlastní produkční fáze (NpT)
22
Analýza MD trajektorií redukce informační exploze
snaha získat snadno uchopitelné a obsažné informace RMSD, Rg, strukturní parametry, RDF difúzní koeficienty essenciální dynamika (PCA) termodynamické veličiny z
sam. přednáška (Tomáš Kubař)
Analýzy – ukázka RMSD
NESTABILNÍ
kvadruplex bez iontů (modře), nativní kvadruplex (červeně) díky Naďo
Analýzy – ukázka B-faktory Bi = 8/3 π2 u2i
23
Esenciální dynamika hledání biologicky relevantních pohybů
vzájemně korelovaná pohyby
Quo vadis, MD?
T. Sprik
Quo vadis? výš – větší a komplexnější systémy (membránové pumpy, komplexy protein/DNA) dál – delší časové škály, pozorovat biologicky relevantní události (folding, recognition ...) a lépe zkoumat fázový prostor (lepší odhady entropických veličin) rychleji – vyšší výkon PC
24
Literatura MacKerell, Jr. A. Empirical Force Fields ... J. Comput. Chem. 25: 1584-1604, 2004
25
