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
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
0.74 = r0 vazebná vzdálenost
( ) (r − r ) + 1 ∂ E (r ) (r − r ) 2!
1 ∂E r0 E (r ) = E ro + 1! ∂r
( )
E (r ) =
(
k r − ro 2
)
2
2
2
0
0
∂r
2
0
+ ...
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
(
k E (r ) = r − ro 2
)
Hookův zákon silová konstanta
)
2
∂ 2 E (r ) ∂r
2
=k
tato aproximace nedovolí disociaci vazeb!
r0
r
(
)
2 1 − a ( r − r0 ) E = ED 1 − e 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
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 H HC N H
H
N
Ala
HC
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
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
⎞ ⎟ ⎟ ⎠
fitováním energetických hyperploch vypočtených referenční QM metodou
meff =
m1 m2 m1 + m2
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°
Deformace torzí
2.9 kcal/mol
E
HH
60
0
H H HH
H
120
240
Degrees of Rotation θ H
H
H H
180
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)
(
(
k i ,t E=∑ 1 + cos ni φi − φi , 0 2 i
))
Deformace torzí
C-C-C-C H-C-C-H
Nepřímé torze např. pro popis vhodné geometrie aminoskupin
AMBER 2-3-1-4 fáze 180° n=2 CHARMM
(
k E = ω ω − ω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 = uatr R i , R j + urep R i , R j
j
k
)
více celá lekce: Slabé interakce (Petr Jurečka)
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
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 polarizaci lze zavést dodatečně
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
Molekulová mechanika (
)
(
)
2 kr Eb = r − r0 2 2 k Ea = θ θ − θ 0 2
(
(
kt 1 + cos nφ − φ0 2 1 qi q j Ec = 4πε 0 ε r rij
Et =
Evdw
⎛ σ ij = −2ε ij ⎜ ⎜r ⎝ ij
6
))
⎛σ ⎞ ⎟ + ε ij ⎜ ij ⎜r ⎟ ⎝ ij ⎠
12
⎞ ⎟ ⎟ ⎠
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
0.000
0.00000
3
DUMM
DU
M
2
1
0
1.522
111.100
4
N
N
M
3
2
1
1.335
116.600
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
10
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
0.000
0.00000
H 12 13
180.000 -0.41570
0.000 -0.56790
torze
H9
7
parc. náboj
8 6
H10
O H11 N H
5
4
Ala
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
snížení počtu nekov. interakcí – zavádí se cutoff
3 4
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)
Cutoff
T. Sprik
Cutoff pro elektrostatiku? elektrostatická interakce dvou monopólů vyhasíná pomalu (r–1)
na vzdálenosti 10r činí 10% vzhledem k interakci na vzdálenosti r
lze řešit trikem při použití periodických okrajových podmínek (PBC) – Ewaldovou sumací
Další zjednodušení popis celých skupin (united atom)
např. skupina –CH3 se bude popisovat jako jeden „pseudoatom“
Běžná silová pole organické molekuly
MM2, MM3, CVFF ...
biomakromolekuly
AMBER z
parm94, 98, 03
CHARMM OPLS+
Molekulová dynamika klasická molekulová dynamika
Molekulová dynamika časový vývoj systému
Molekulová dynamika ∂ ri 2
Fi = mi a i = mi
∂t
2
= mi &r&i
integrace Newtonových rovnic
∂E ri Fi = − ∂ri ri potenciál MM
1 2 Fi (t ) stará ri (t + Δt ) = ri (t ) + (Δt )v i (t ) + (Δt ) 2 mi nová stará stará
aktualizované souřadnice v čase Δt (1-2 fs)
velocity Verlet
Aktualizace rychlostí Fi (t ) 1 ⎛ Δt ⎞ v i ⎜ t + ⎟ = v i (t ) + (Δt ) 2 mi 2⎠ ⎝
rychlosti se aktualizují ve dvou krocích Δt/2
Fi (t + Δt ) ⎛ Δt ⎞ 1 v i (t + Δt ) = v i ⎜ t + ⎟ + (Δt ) 2⎠ 2 mi ⎝
teplota a primitivní realizace termostatu N N 1 1 3 2 2 mi v i ≅ Nk BT T= mi v i ∑ ∑ 2 3 Nk B i =1 i =1 2
Trajektorie
stav 0 ri0
stav 1 ri1
stav 2
stav 3
Fi = mi ai
snímky = MD trajektorie (ri,j, vi,j) termostat (T)
díky Evo
Jak začít? generace náhodných rychlostí pro danou teplotu
⎛ v i2 ⎞ 1 p( v i ) = exp⎜ − 2 ⎟, N 0, σ 2 ⎜ 2σ ⎟ σ 2π ⎝ ⎠ k BT 2 2 2 2 σ = vi − vi = vi = mi
(
)
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 více přednáška (Dan Svozil)
periodické okrajové podmínky
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 ...
Solut - protein
Solut + protiionty
Solut + protiionty + voda (box)
93.678 atomů
Periodické okrajové podmínky
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)
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
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
